Maximum-principle preserving space-time isogeometric analysis

In this work we propose a nonlinear stabilization technique for convection-diffusion-reaction and pure transport problems discretized with space-time isogeometric analysis. The stabilization is based on a graph-theoretic artificial diffusion operator and a novel shock detector for isogeometric analysis. Stabilization in time and space directions are performed similarly, which allow us to use high-order discretizations in time without any CFL-like condition. The method is proven to yield solutions that satisfy the discrete maximum principle (DMP) unconditionally for arbitrary order. In addition, the stabilization is linearity preserving in a space-time sense. Moreover, the scheme is proven to be Lipschitz continuous ensuring that the nonlinear problem is well-posed. Solving large problems using a space-time discretization can become highly costly. Therefore, we also propose a partitioned space-time scheme that allow us to select the length of every time slab, and solve sequentially for every subdomain. As a result, the computational cost is reduced while the stability and convergence properties of the scheme remain unaltered. In addition, we propose a twice differentiable version of the stabilization scheme, which enjoys the same stability properties while the nonlinear convergence is significantly improved. Finally, the proposed schemes are assessed with numerical experiments. In particular, we considered steady and transient pure convection and convection-diffusion problems in one and two dimensions.



page 15

page 17


On differentiable local bounds preserving stabilization for Euler equations

This work is focused on the design of nonlinear stabilization techniques...

An experimental comparison of a space-time multigrid method with PFASST for a reaction-diffusion problem

We consider two parallel-in-time approaches applied to a (reaction) diff...

An Eulerian-Lagrangian Runge-Kutta finite volume (EL-RK-FV) method for solving convection and convection-diffusion equations

We propose a new Eulerian-Lagrangian Runge-Kutta finite volume method fo...

Maximum Principle Preserving Finite Difference Scheme for 1-D Nonlocal-to-Local Diffusion Problems

In a recent paper (see [7]), a quasi-nonlocal coupling method was introd...

Convergence analysis of the variational operator splitting scheme for a reaction-diffusion system with detailed balance

We present a detailed convergence analysis for an operator splitting sch...

Maximum bound principles for a class of semilinear parabolic equations and exponential time differencing schemes

The ubiquity of semilinear parabolic equations has been illustrated in t...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

Many different applications in science and industry require solving problems satisfying some sort of positivity or maximum principle (MP) property. These include scalar transport problems, compressible flows, or fluid-based MHD simulations, among others. These problems are of particular interest in a variety of industries and scientific research areas, such as the chemical industry, aviation, aerospace, or nuclear fusion research, just to cite few examples.

Some of these problems exhibit a multiscale nature in time. In those cases, explicit methods are not suitable, since the smallest time scales pose very stringent stability conditions to the time step length, i.e., fully resolved time simulations are required. Thus, implicit methods are favored in applications where the smallest time scales are not of scientific or engineering interest.

As a result, schemes that preserve monotonicity (or at least positivity) for implicit time integration are of special interest. The standard technique to attain such schemes is adding nonlinear artificial diffusion (usually called shock capturing). The common ingredients of a shock capturing or nonlinear stabilization method are the following. The first ingredient is the artificial diffusion, which needs to be sufficient to eliminate non-physical oscillations. The schemes in [9, 8, 4, 5] use an element-based artificial diffusion with a standard PDE-based diffusion operator. The drawback of this choice is the fact that the discrete maximum principle (DMP) only holds under unpractical mesh restrictions. This problem has been solved by Guermond and Nazarov in [15, 12] by replacing the PDE-based diffusion operator by an edge or graph-theoretic diffusion operator; see [2, 3, 19, 22, 23] for schemes that preserve the DMP on arbitrary meshes using a graph-Laplacian. The second ingredient is a shock detector, which is the term responsible of deactivating the artificial diffusion in smooth regions. A good shock detector is of vital importance for minimizing the numerical diffusion while satisfying a DMP. One example of shock detector is the one developed in 1D by Burman in [8] and later extended to multiple dimensions by Badia and Hierro [4]. The last ingredient consists on perturbing the mass matrix. One option is a full lumping of the mass matrix, but it can lead to unacceptable phase errors. Instead, a nonlinear lumping is used, e.g., in [3, 2], using the same shock detector to lump the mass matrix. Other alternatives can be found in [19, 14]. It is worth mentioning that all previous stabilization methods yield a very stiff nonlinear system of equations. In fact, some of the methods proposed in the literature are not even Lipschitz continuous and thus ill-posed (see [7]). In practice, the nonlinear convergence of these methods is unacceptably slow, making hard its practical use. To solve this problem, Badia et al. [3, 2] have designed differentiable nonlinear stabilization terms, noticeably improving the nonlinear convergence.

The methods commented above have an algebraic nature and provide some type of DMP for the nodal values. The monotonicity of the nodal values only translates into monotonic solutions if the FE space satisfies the convex hull property, which is only true in the first order case. As a result, using the ideas above it does not seem possible to design monotonic second or higher order methods. Recently, Kuzmin and coworkers [1, 22], have proposed instead the usage of Bernstein–Bèzier finite elements, since they satisfy the convex hull for high-order. However, the temporal dimension is discretized using Backward-Euler or strong stability preserving (SSP) Runge-Kuta methods (see [16]). In the first case, the problem is first order in time, whereas in the second case, a CFL-like condition arises [21], since high-order SSP methods pose a restriction on the time step size similar to the ones in explicit methods [16].

The main contribution of the present work is the development of a high-order (both in space and time) and DMP-preserving discretization for the convection-diffusion-reaction and pure transport problems. This is achieved by combining the nonlinear stabilization techniques in [3, 2] with a new shock detector for arbitrary order space-time isogeometric analysis. Another novelty of this study is the stabilization in the time direction, which is performed in a similar manner as in space. This results in an unconditionally stable high-order method in time (and space). However, the space-time method requires to solve the whole space-time problem at once, which increases the computational cost. Hence, we also propose a partitioned approach in the temporal direction, where one can determine the width of the time slab to be computed every time. This strategy allows us to maintain a reasonable computational cost while having a high-order scheme in space and time, as well as satisfying the DMP without any CFL-like condition. Finally, we also propose a differentiable version of the above scheme. This allow us to use Newton’s method, which improves nonlinear convergence significantly. The method proposed in the present work has been implemented and tested making use of the FEMPAR library [6].

This paper is structured as follows. First, we introduce the problem, its discretization, and monotonicity properties for scalar problems in Sect. 2. Then, the stabilization techniques are introduced in Sect. 3. Sect. 4 is devoted to the partitioned time integration scheme. Afterwards, we introduce a regularized version of the stabilization term in Sect. 5. Finally, we show numerical experiments in Sect. 6 and draw some concluding remarks in Sect. 7.

2. Preliminaries

2.1. Convection-Diffusion problem

We consider a transient convection-diffusion problem with Dirichlet boundary conditions. Let be a -cube, where is the number of spatial dimensions. Then, the problem reads:


where is a divergence-free convection velocity, is a scalar constant diffusion, and is the body force. In the case of pure convection (), boundary conditions are only imposed at the inflow , where

is a unit vector outward-pointing normal to the boundary. We also define the outflow boundary as

. Moreover, we will also consider the steady problem, which is obtained by dropping the time derivative term and the initial conditions. It is important to mention that a reaction term can be included without harming any of the properties satisfied by the schemes introduced below. However, a convection–diffusion–reaction problem only satisfies a MP if the minimum is negative and the maximum positive (analogously for its proposed discretizations). In other words, it only satisfy a weak MP, see [7]. In order to simplify the discussion below, we will limit the present work to pure convection and convection–diffusion problems.

2.2. Discretization

In this work, we consider a standard B-spline discretization with interpolative boundaries (see

[10]). A spline of order in the variable is a piecewise polynomial function in of degree . The values of in which different polynomials meet are call knots. Knots might be placed at the same location, i.e. can be repeated. When the knots are not repeated, the first derivatives of the spline are continuous. When a knot is repeated times, only the first derivatives are continuous across that knot. Knots are sorted in increasing order and collected in the so called knot vector . Given a knot vector, B-splines of order are basis functions for spline functions of the same order. B-splines are constructed in a recursive way using the Cox-de Boor formula:


for . By construction, is compact, non-negative, and non-zero in .

Let us consider the domain and the uniform partition into sub-intervals of size . The open knot vector is defined as follows. The first knots are located at zero, i.e., . The last knots are located at , i.e., . The interior points are equi-distributed, with , for . It leads to a basis (for ) for a space of splines in and a partition of unity, i.e., for . Any spline of order in can uniquely be defined by the control points as the linear combination of B-splines . In one dimension, the basis functions obtained from an open knot vector are interpolatory at the extremes, i.e., and (see Fig. 1). For a first order polynomial in , it holds , where are called the Greville abscissae [11, 24].

Let us consider the number of partitions per dimension with , for . We represent with the set of multi-indices with . Every can be expressed as , where is the spatial index and is the temporal index. The

-dimensional B-spline is defined as the tensor product of

unidimensional B-splines . Notice that a Greville abscissa in the case of a multidimensional spline reads .

We define the space of splines . We use the notation . The order is omitted since it is assumed to be fixed. Thus, every spline can be written as . Furthermore, we define the following sets of indices, which are useful for the definition of the forthcoming schemes. The set of neighbors of is defined as . We define as the set of indices whose associated shape functions intersect with the support of .

Figure 1. Representation of the basis functions of in one dimension, with its associated Greville abscissae.

We use standard notation for Sobolev spaces. The scalar product is denoted by for . However, we omit the subscript for . The norm is denoted by .

2.3. Discrete problem

The weak form of (1) using the Galerkin method reads: find such that on , on , and


where and are projections of and to , respectively, such that the local DMP is satisfied (see Def. 2.2). Furthermore, we can rewrite the previous discrete problem in matrix form as , where , and for . Notice that we have not applied the boundary conditions yet. To apply boundary conditions the space of test functions is restricted to , and the force vector is redefined as .

2.4. Monotonicity properties

In this section we define all the properties that we demand our scheme to fulfill. In this case, since we are using a space-time discretization, it becomes more useful to define these properties in a space-time sense. This means that the variation of in the temporal direction will also be taken into account to define an extremum. Hence, we define the concept of a local discrete extremum as follows.

Definition 2.1 (Local Discrete Extremum).

The function has a local discrete minimum (resp. maximum) on , if (resp. ) .

It is also important to define the concepts of local and global space-time DMP. The latter is a slightly weaker property than the former, but it is more useful in the present study. A local DMP is a stronger property because it implies that no oscillations can appear, while the global DMP only implies that the global extrema are located at the boundary conditions.

Definition 2.2 (Local space-time Dmp).

A solution satisfies the local discrete maximum principle if for every

Definition 2.3 (Global space-time Dmp).

A solution satisfies the global discrete maximum principle if the global extrema are located at boundary conditions, i.e., for every


Finally, let us recall the definition of linearity-preservation, which is a desired property to achieve high-order convergence in smooth regions (see [20]).

Definition 2.4 (Linearity-preservation).

A stabilization term, , is said to be linearity-preserving if, for a solution that is linear in any direction, then the stabilization term becomes null, i.e., for .

3. Lipschitz-continuous nonlinear stabilization

In this section we define a nonlinear stabilization operator, , to be added to the discrete problem (3), such that it satisfies at least the global DMP in Def. 2.3. Let us define . We also enforce that, for any ,

  1. has compact support: if ,

  2. is symmetric: ,

  3. is conservative: .

In order to achieve these requirements, we recall the stabilization term in [2], which is defined as


for any . Here, is the graph-Laplacian operator defined as (see [13, 2]), and is the artificial diffusion defined as


We denote by the shock detector used for computing the artificial diffusion parameter. The idea behind the definition of this detector is to ensure that the global DMP defined in Def. 2.3 is satisfied using a minimal amount of artificial diffusion, i.e., the lower admissible value of . A shock detector must be a positive real number, which takes value 1 when is an inadmissible value of (i.e., local discrete extremum) and smaller than 1 otherwise; to have linearity preservation (see Def. 2.4), it must be equal to 0 for . In this section, we propose an isotropic approach for , which consists in using the shock detector in [2] in all directions (including time).

Figure 2. Representation of the polytope in two dimensions, the symmetric node of with respect to , and .

In order to introduce the shock detector, let us recall some useful notation from [2]. Let be the vector pointing from Greville abscissae to with and . Let us take the set of Greville abscissae for as vertices of a polytope in dimensions. In particular, let us name this polytope . Let be the point at the intersection between and the line that passes through and that is not (see Fig. 2). The set of all for all is represented with . We define . Given in two dimensions, let as call and the indices of the vertices such that they define the edge in that contains . We define as the linear interpolation of and at , i.e. . For higher dimensions, is defined analogously. Given the facet of where lies, is the linear interpolation at of the control points whose Greville abscissae are at the same facet.

Notice that it is essential to use Greville abscissae since they satisfy that for a linear function , . Therefore, one can construct easily linear approximations of the unknown gradients that are exact for . Furthermore, one can define the jump and the mean of a linear approximation of the unknown gradient at Greville abscissa in direction as


In the present work we will use the same shock detector developed in [2], which reads


From Lm. 3.1 in [2] we know that (11) is valued between 0 and 1, and it is only equal to one if is a local discrete extremum (in a space-time sense as in Def. 2.1). Since the linear approximations of the unknown gradients are exact for , the shock detector vanishes when the solution is linear in all dimensions (including time). This result follows directly from [2, Th. 4.5].

Supplementing the discrete problem (3) with the above stabilization term, the stabilized problem reads: Find such that on , at , and


which in turn can be expressed in matrix form as


where for .

Theorem 3.1 (Dmp).

The solution of the discrete problem (13) using the shock detector (11) satisfies the global DMP in Def. 2.3 if and, for every control point such that is a local discrete extremum, it holds:


Moreover, the resulting scheme is linearity-preserving as defined in Def. 2.4, i.e. for .


Let us assume that is a discrete maximum. Then, (13), for and before applying boundary conditions, reads


Therefore, can be computed as


Since is an extremum, which implies , the stabilization term ensures (14) by construction. Hence, the coefficients that multiply are in , and the sum of all these coefficients add up to one. Therefore, is a convex combination of its neighbors (including boundary conditions ). Since is a maximum and a convex combination of its neighbors, then for some . Furthermore, it can also be proved that is a convex combination of all its neighbors but , and vice versa is a convex combination of all its neighbors but . Hence, by induction, we know that extrema at any control point are bounded by the boundary conditions. Thus, the global DMP is satisfied.

From [2, Th. 4.5], implies for any . By definition, the stabilization term also vanishes if (see (6) and (7)). Therefore, the scheme is linearity-preserving as defined in Def. 2.4. ∎

Theorem 3.2.

The diffusion defined in (7) introduces the minimal amount of artificial dissipation such that condition (14) is satisfied when .


The proof follows the same lines as in [2, Th. 4.4]. We do not include it for the sake of conciseness. ∎

Finally, we proof Lipschitz continuity of the stabilization term. In this particular case, the proof follows the same reasoning as in [2]. We do not include all details for the sake of conciseness and refer the reader to the previously cited work.

Theorem 3.3.

Let be a non-degenerate partition, then with the shock detector (11) is Lipschitz continuous in for and a bounded , since


is satisfied.


The proof is analogous to the one in [2, Th. 6.1]. The only difference arises from the bound for


In this case, using Cauchy-Schwarz inequality, the inverse inequality for , and , we get


where is the number of spatial dimensions, is the distance between knots for the spatial directions and is the distance between knots for the time direction. ∎

4. Time partitioned scheme

Hitherto, we have only considered the solution of the whole space-time problem at once. In order to substantially reduce the computational cost, we propose the division of the time integration in several time subdomains, considering the proposed space-time formulation at every subdomain. Namely, the problem set in , will be decomposed in for , with and . We define the length of the subdomain as , and restrict its possible values to for some , where is the order of the spline space, and is the distance between knots in the temporal direction. Notice that we are only using discretizations formed from the tensor product of discretizations in 1D. Therefore, with the particular choice of , will always be the temporal coordinate of a layer of knots. Hence, performing this kind of partitions is straightforward. Other partitions might be considered, however we choose the previous one because it is particularly simple to use it in our implementation.

The approximation space of splines for every subdomain is obtained as follows. Given the complete domain , we discretize it as described in Sect. 2.2, resulting in a spline space . Then in order to reduce the coupling between partitions, we insert knots at . The resulting spaces at each subdomain, say , are fully decoupled. However, due to causality in time there exists a sequential coupling between subdomains, i.e. the information travels in the positive direction. In other words, the solution at subdomain will affect the solution at , but not the opposite. Therefore, we impose that the initial conditions at subdomain are equal to the solution at the final time of subdomain , i.e. . After imposing this restriction, the complete approximation space, , is , and coupled sequentially. Hence, each subdomain can be solved sequentially, and thus the computational cost is significantly reduced.

The partitioned space-time scheme with nonlinear stabilization reads as follows. For ; find such that on , with , and


Due to the partition will only be piecewise continuous in time. Let us proof now that the scheme still satisfies the global DMP.

Lemma 4.1.

The solution of problem (21), using the shock detector defined in (11), satisfies the global DMP (Def. 2.3) if in and, for every control point such that is a local discrete extremum, conditions (14) hold.


From Th. 3.1 it is easy to see that conditions (14) hold for the first subdomain. Hence, the solution at the first subdomain at time , , is bounded by the initial and boundary conditions. Since are the initial conditions for the second subdomain, then again from Th. 3.1 it is known that the solution in the second subdomain is bounded by , and thus by the initial and boundary conditions. Therefore, by induction, we conclude that the global DMP is satisfied in the whole domain. ∎

5. Differentiable stabilization

In this section, we introduce a different version of the previous operators. As exposed in [3, 2], the regularization of all non-differentiable operators in the stabilization term improves the nonlinear convergence, and allows us to use Newton’s method. We use the same strategy introduced in [2]. Then, the shock detector reads


where is a parameter to prevent division by zero, and the regularized absolute values by


Notice that . With this regularization, the quotient in the shock detector might become greater than one, thus we need to smoothly limit its value to one. To this end we recall , which reads


and clearly is twice differentiable and bounded above by 1. The differentiable version still satisfies the requirements for a shock detector, i.e., it is a real value in and it is equal to 1 if is a local extrema. This result follows directly from [2, Lm 7.1].

Furthermore, for the definition of the artificial diffusion we need to regularize the maximum function. We choose again the same strategy as in [2], and define as


Finally, we can define the twice differentiable artificial diffusion parameter as


and the stabilization operator reads


In order to obtain a differentiable operator, we have added a set of regularizations that rely on different parameters, e.g., . Giving a proper scaling of these parameters is essential to recover theoretic convergence rates. In particular, we use the following relations


where is the spatial dimension of the problem, is a characteristic length, and and are of the order of the unknown.

6. Numerical experiments

In this section we present some numerical experiments showing the behavior of the scheme previously introduced. First, a convergence analysis is performed in order to assess the correctness of the proposed scheme and its implementation. Then, we assess the performance of the proposed stabilization method for high-order discretizations, including a brief analysis of the effect of the regularization.

6.1. 1D Transient Diffusion

The purpose of this test is assessing the partitioned time integration scheme in Sect. 4. To this end, we solve the following problem for and ,


where . This problem has as exact solution. We perform a convergence analysis where the mesh is successively refined in the time direction for first, second, and third order discretizations. In particular, the distance between knots in the temporal direction is refined as for the first order discretization. In spatial directions, the distance is small enough to prevent that spatial discretization errors affect the analysis. Second and third order discretizations are obtained using the following -refinement (see [10] for more details). We refine the discretization such that the number of control points increase at the same rate as a Lagrangian FE discretization does when its order is increased. Fig. 3 shows the result of -refinements to and discretizations, for an interior subset of the discretization. Henceforth, we will use this kind of -refinement in order to increase the discretization order.

Figure 3. Second and third order discretizations obtained from the -refinement of an initial first order discretization. Notice that shape functions are depicted for interior knots, at boundary knots shape functions become interpolatory, see Fig. 1.
(a) Relative norm of the error.
(b) Relative norm of the error.
Figure 4. Convergence in time results for problem (29), using standard and partitioned space-time schemes.

We measure the relative error in norm and semi-norm, and compute the resulting convergence rate. Errors in norm and semi-norm are depicted in Fig. 4 (a) and (b), respectively. In Table 1 the measured convergence rates are shown for the original non-partitioned scheme and the proposed in Sect. 4. We observe a slight increase in the error for the partitioned scheme. However, the obtained results show optimal convergence rates for both schemes introduced above.

Order Method convergence convergence
1 p-ST -1.98 -0.98
1 ST -2.04 -0.96
2 p-ST -3.03 -2.00
2 ST -3.00 -2.01
3 p-ST -3.99 -3.00
3 ST -3.99 -2.99

p-ST: Partitioned space-time, ST: space-time.

Table 1. Measured convergence rates in norm and semi-norm, for problem (29).

6.2. Steady convection

In this experiment we assess the convergence of the stabilized schemes introduced in Sect. 5. We use a steady pure convection problem with a non-monotonic smooth solution. In particular, we solve the following problem for ,


where , , and . The analytical solution of the above problem reads . The convergence analysis is performed for first, second, and third order discretizations. We use a standard nonlinear solver (see [3] for details), with a nonlinear tolerance . The following stabilization parameters have been used: , , , . The selection of these parameters is based on the outcome of previous works [3, 2] and Sect. 6.3.

Convergence plots are shown in Fig. 5 and the corresponding convergence rates in Table 2. As expected, it is observed that the scheme recovers second order convergence in the error norm and first in the error semi-norm. It is known that the stabilized scheme should recover second order convergence for . However, due to peak clipping errors, higher convergence rates are not expected even if a higher order discretization is used [18]. In any case, we do observe that the error diminishes as the discretization order is increased using the -refinement previously defined.

(a) Relative norm of the error.
(b) Relative semi-norm of the error.
Figure 5. Convergence in space results for problem (30).
Order convergence convergence
1 1.77 0.88
2 1.85 0.96
3 1.86 0.99
Table 2. Measured convergence rates in and norms, for problem (30).

6.3. Nonlinear convergence

In the current test, we aim to briefly analyze the effect of the stabilization parameters on the nonlinear convergence of the method. To this end, we solve the following 1D pure convection problem with discontinuous initial conditions.


where , , , and , where is the well-known zero-centered Heaviside function. First and second discretization orders are used in a coarse mesh of control points. To obtain the second order mesh, we perform the -refinement as in the previous experiment.

We refer the reader to [3, 2] for a deeper analysis on the effect of each regularization parameter. Therein, the same family of shock detectors is used in the context of first order cG and dG Lagrangian FEs. In the present study, we analyze the effect of the regularization globally using a fixed relation between the different parameters. In particular, we use the following parameters: , , , where . Furthermore, the effect is also compared as is incremented, particularly for . In addition, the non-regularized version is also used to show the improvement in the nonlinear convergence. The relaxed Picard and hybrid nonlinear solvers presented in [3] are used, the nonlinear tolerance is set to , and a maximum of 300 nonlinear iterations are allowed.

Fig. 6 and 7 show the results for first and second order discretizations, respectively. In general terms, as is increased or is decreased, sharper solutions are observed. However, nonlinear iterations increase. As expected, the hybrid method outperforms the relaxed Picard method. Even though it requires more nonlinear iterations, the non-regularized detector might be a simpler (parameter-free) alternative to the regularized one. Finally, it is worth mentioning that a slight increase in the required number of iterations is observed as the discretization order is increased. However, the obtained results are more accurate, i.e., the discontinuity becomes sharper.

Figure 6. Effect of the regularization parameters for first order discretizations. The numbers in legends are the number of nonlinear iterations performed. First number is for relaxed Picard and the next for hybrid scheme, both for the regularized stabilization. The number in brackets is the number of iterations required to converge the non-differentiable method using relaxed Picard scheme.
Figure 7. Effect of the regularization parameters for second order discretization. The numbers in legends are the number of nonlinear iterations performed. First number is for relaxed Picard and the next for hybrid scheme, both for the regularized stabilization. The number in brackets is the number of iterations required to converge the non-differentiable method using relaxed Picard scheme.

6.4. 1D Sharp layer propagation

The performance of the stabilization schemes is analyzed as the discretization order is increased. To this end, we use again the previous problem (31). The regularization parameters are kept fixed, while the discretization is modified both in terms of the order of accuracy and the number of control points.

We use a nonlinear tolerance of , and limit to 300 the maximum nonlinear iterations. The regularization parameters used are , , , and . With this setting, we solve the above problem using a discretization that keeps the number of control points fixed as the order is increased, and another one using the -refinement defined in the previous experiment. For the former, we use a discretization of 120 by 60 control points. For the latter, we start with a first order discretization of 120 by 60 control points and refine as previously explained.

In Fig. 7(a) the solution at is shown for different orders and fixed number of control points, and using the -refinement in Fig. 7(b). We observe that for non-smooth solutions, fixing the number of control points and increasing the order does not improve the results. In the case of Fig. 7(b), as expected, we observe better approximations as the order is increased using the -refinement. Hence, better results might be expected as the order is increased for problems that combine discontinuities and smooth profiles.

In Fig. 9, similar results are shown when using the time integration scheme proposed in Sect. 4. A small degradation of the results can be seen in Fig. 8(a) as we increase the discretization order. In a similar trend, we observe less improvement in Fig. 8(b) than in Fig. 7(b). We attribute this degradation to the time partitions, which becomes more evident as the subdomains are smaller.

(a) Solutions increasing the order while keeping fixed the number of control points.
(b) Solution increasing the order using the -refinement process described above.
Figure 8. Solution of problem (31) at for first to fourth order discretizations.
(a) Solution using fixed number of control points, and time integration defined in Sect. 4 with 5 partitions.
(b) Solution using -refinement, and time integration defined in Sect. 4 with 5 partitions.
Figure 9. Solution of problem (31) at for first to fourth order discretizations.

6.5. Boundary layer

In this section the effect of the discretization order in a convection-diffusion problem is analyzed. To this end, we solve a problem with the propagation of a sharp layer and a boundary layer. In particular, we solve for


where , , and the boundary conditions are defined as


For this test, we use the following settings: a nonlinear tolerance of , , , , and . In Fig. 9(a), the solution for is depicted. The converged solution does not exhibit any oscillation. Very sharp layers are obtained for this parameter setting. In Fig. 9(b), we show the profile of the solution at for different orders. In this case, we start with a discretization of 50 control points per direction. Then, we increase the order using the -refinement used previously. As previously observed for transient problems, we observe an improvement of the solution as the order is increased.

(a) 3D representation of the solution.
(b) Profiles at .
Figure 10. Solution of problem (32) using scheme (12), and different discretization orders.

6.6. Three Body rotation

Finally, we solve the transient pure convection problem (31) in for , with . Initial conditions are given in [17]. Its interpolation in a first order control point mesh is depicted in Fig. 10(a). The analytical solution of this problem is simply the translation of the profiles in the direction of the convection. In particular, for , one revolution is completed and the solution is equal to the initial conditions. The purpose of this test is to evaluate how diffusive is the proposed scheme. We perform this evaluation evolving the solution until and comparing the results with the initial conditions.

The solution is computed using scheme (12) in combination with the shock detector in (22). We use the following parameters for the stabilization: , , , and . Different meshes, time partitions, and discretization orders are used in this experiment. We start with a linear discretization of control points in space, and 500 in time divided in 125 subdomains. Then, we increase the discretization order to using the -refinement. In order to compare first and second order discretizations, but using a similar number of control points we use a discretization with control points in space, and 1000 in time divided in 250 subdomains. Finally, we assess the effect of the partitions in the temporal direction. We compare the previous discretization of control points divided in 125 subdomains, with the same discretization divided in 250 subdomains. We do the same comparison for the second order discretization using 125 subdomains and when it is divided in 250 subdomains.

(a) Initial conditions of the 3 body rotation.
(b) Profile at .
Figure 11. Three body rotation test initial conditions.

Fig. 13 shows the solutions for meshes, and 125 subdomains in time, whereas Fig. 14 show the ones for 250 subdomains. In both cases, a great improvement can be observed as we increase the discretization order. However, the computational cost is also increased. It is interesting to compare the solutions for first and second order discretizations using meshes with similar amount of control points, namely solutions at Fig. 12 and 12(b). For this particular problem, using a higher order discretization with similar number of control points does not improve the solution, which it is actually slightly more diffusive for . It is also worth mentioning that increasing the discretization order does not modify the behavior of the solution in terms of clipping or terracing. Comparing Figs. 13(a) and 12(a), we observe that the scheme becomes more dissipative as the number of partitions is increased. This is even clearer in Fig. 15, where the profile of the solution at is depicted.

(a) 3D view.
(b) Profile at .
Figure 12. Three body rotation test results at using scheme (12), , , , and . A first order discretization of control points is used with 250 subdomains in the temporal direction.
(a) Solution for a first order discretization.
(b) Solution after one -refinement.
Figure 13. Three body rotation test results at using scheme (12), , , , and . A first order discretization of control points is used. The second order discretization is obtained using -refinement. 125 subdomains in the temporal direction have been used.
(a) Solution for first order discretization.
(b) Solution after one -refinement.
Figure 14. Three body rotation test results at using scheme (12), , , , and . A first order discretization of control points is used. The second order discretization is obtained using -refinement. 250 subdomains in the temporal direction have been used.
(a) Solution for first order discretization.
(b) Solution after one -refinement.
Figure 15. Three body rotation test profiles for at using scheme (12), , , , and . A first order discretization of control points is used. The second order discretization is obtained using -refinement. 125 and 250 subdomains in the temporal direction have been used.

7. Conclusions

In the present work, an extension of [2] to isogeometric analysis methods have been developed. The proposed method is unconditionally DMP preserving for arbitrary high-order discretizations in space and time without any CFL-like condition. Furthermore, it is shown to be linearity-preserving in a space-time sense. Moreover, the regularized version is shown to yield better convergence behavior, specially when for the hybrid Picard-Newton method.

Moreover, the numerical experiments show that increasing the discretization order yield much better solutions. However, as the order is increased the number of control points and the computational cost is also increased. On the contrary, if the order is increased while the number of control points is fixed, then similar or even slightly more diffusive results are obtained for non-smooth solutions. Hence, for problems with regions of smooth and non-smooth solutions a high-order is expected to outperform linear discretizations with similar amount of control points. Furthermore, a method capable of providing solutions that satisfy the DMP for high-order discretizations is of special interest in -adaptive schemes, since the usage of first order discretizations in shocks is not required.

In addition, a partitioned scheme that does not harm any monotonicity property is presented. This scheme reduces significantly the computational cost of the original space-time scheme. It is important to mention, that this partitioning slightly increases the error. However, this approach allows finer meshes, and thus, in practice better solutions can be obtained.


J. Bonilla gratefully acknowledges the support received from ”la Caixa” Foundation through its PhD scholarship program. S. Badia gratefully acknowledges the support received from the Catalan Government through the ICREA Acadèmia Research Program. We acknowledge the financial support to CIMNE via the CERCA Programme / Generalitat de Catalunya.


  • [1] R. Anderson, V. Dobrev, T. Kolev, D. Kuzmin, M. Quezada de Luna, R. Rieben, and V. Tomov, High-order local maximum principle preserving (MPP) discontinuous Galerkin finite element method for the transport equation, Journal of Computational Physics, 334 (2017), pp. 102–124.
  • [2] S. Badia and J. Bonilla, Monotonicity-preserving finite element schemes based on differentiable nonlinear stabilization, Computer Methods in Applied Mechanics and Engineering, 313 (2017), pp. 133–158.
  • [3] S. Badia, J. Bonilla, and A. Hierro, Differentiable monotonicity-preserving schemes for discontinuous Galerkin methods on arbitrary meshes, Computer Methods in Applied Mechanics and Engineering, 320 (2017), pp. 582–605.
  • [4] S. Badia and A. Hierro, On Monotonicity-Preserving Stabilized Finite Element Approximations of Transport Problems, SIAM Journal on Scientific Computing, 36 (2014), pp. A2673—-A2697.
  • [5] S. Badia and A. Hierro, On discrete maximum principles for discontinuous Galerkin methods, Computer Methods in Applied Mechanics and Engineering, 286 (2015), pp. 107–122.
  • [6] S. Badia, A. F. Martín, and J. Principe, FEMPAR: An Object-Oriented Parallel Finite Element Framework, Archives of Computational Methods in Engineering, 25 (2018), pp. 195–271.
  • [7] G. R. Barrenechea, E. Burman, and F. Karakatsani, Edge-based nonlinear diffusion for finite element approximations of convection–diffusion equations and its relation to algebraic flux-correction schemes, Numerische Mathematik, (2016), pp. 1–25.
  • [8] E. Burman, On nonlinear artificial viscosity, discrete maximum principle and hyperbolic conservation laws, BIT Numerical Mathematics, 47 (2007), pp. 715–733.
  • [9] E. Burman and A. Ern, Nonlinear diffusion and discrete maximum principle for stabilized Galerkin approximations of the convection–diffusion-reaction equation, Computer Methods in Applied Mechanics and Engineering, 191 (2002), pp. 3833–3855.
  • [10] J. A. Cottrell, T. J. R. Hughes, and Y. Bazilevs, Isogeometric analysis: toward integration of CAD and FEA, Wiley, 2009.
  • [11] C. De Boor, B(asic)-Spline Basics, tech. rep., Madison mathematics research center, Winsconsin University., 1986.
  • [12] J.-L. Guermond and M. Nazarov, A maximum-principle preserving finite element method for scalar conservation equations, Computer Methods in Applied Mechanics and Engineering, 272 (2014), pp. 198–213.
  • [13] J.-L. Guermond, M. Nazarov, B. Popov, and Y. Yang, A Second-Order Maximum Principle Preserving Lagrange Finite Element Technique for Nonlinear Scalar Conservation Equations, SIAM Journal on Numerical Analysis, 52 (2014), pp. 2163–2182.
  • [14] J.-L. Guermond and R. Pasquetti, A correction technique for the dispersive effects of mass lumping for transport problems, Computer Methods in Applied Mechanics and Engineering, 253 (2013), pp. 186–198.
  • [15] J.-L. Guermond, R. Pasquetti, and B. Popov, Entropy viscosity method for nonlinear conservation laws, Journal of Computational Physics, 230 (2011), pp. 4248–4267.
  • [16] D. I. Ketcheson, C. B. Macdonald, and S. Gottlieb, Optimal implicit strong stability preserving Runge–Kutta methods, Applied Numerical Mathematics, 59 (2009), pp. 373–392.
  • [17] D. Kuzmin, A Guide to Numerical Methods for Transport Equations, Friedrich-Alexander-Universität, 2010.
  • [18] D. Kuzmin, S. Basting, and J. N. Shadid, Linearity-preserving monotone local projection stabilization schemes for continuous finite elements, Computer Methods in Applied Mechanics and Engineering, 322 (2017), pp. 23–41.
  • [19] D. Kuzmin and J. N. Shadid, Gradient-based nodal limiters for artificial diffusion operators in finite element schemes for transport equations, International Journal for Numerical Methods in Fluids, (2016), pp. 675–695.
  • [20] D. Kuzmin, M. J. Shashkov, and D. Svyatskiy, A constrained finite element method satisfying the discrete maximum principle for anisotropic diffusion problems, Journal of Computational Physics, 228 (2009), pp. 3448–3463.
  • [21] D. Kuzmin and S. Turek, Flux Correction Tools for Finite Elements, Journal of Computational Physics, 175 (2002), pp. 525–558.
  • [22] C. Lohmann, D. Kuzmin, J. N. Shadid, and S. Mabuza, Flux-corrected transport algorithms for continuous Galerkin methods based on high order Bernstein finite elements, Journal of Computational Physics, 344 (2017), pp. 151–186.
  • [23] S. Mabuza, J. N. Shadid, and D. Kuzmin, Local bounds preserving stabilization for continuous Galerkin discretization of hyperbolic systems, Journal of Computational Physics, 361 (2018), pp. 82–110.
  • [24] M. J. Marsden, An identity for spline functions with applications to variation-diminishing spline approximation, Journal of Approximation Theory, 3 (1970), pp. 7–49.