Maximum-principle preserving space-time isogeometric analysis

TL;DR

This paper introduces a nonlinear stabilization method for space-time isogeometric analysis that preserves the maximum principle, allows high-order discretizations without CFL constraints, and reduces computational costs through a partitioned scheme.

Contribution

It presents a novel graph-theoretic artificial diffusion stabilization technique with a shock detector, unconditionally preserving the discrete maximum principle for arbitrary order.

Findings

The method ensures solutions satisfy the discrete maximum principle unconditionally.

A partitioned space-time scheme reduces computational costs while maintaining stability.

A twice differentiable stabilization scheme improves nonlinear convergence.

Abstract

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

Tables

Table 1. Table 1 . Measured convergence rates in L 2 superscript 𝐿 2 L^{2} norm and H 1 superscript 𝐻 1 H^{1} semi-norm, for problem ( 33 ).

Order	Method	$L^2$ convergence	$H^1$ 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

Table 2. Table 2 . Measured convergence rates in L 2 superscript 𝐿 2 L^{2} and H 1 superscript 𝐻 1 H^{1} norms, for problem ( 34 ).

Order	$L^2$ convergence	$H^1$ convergence
1	1.77	0.88
2	1.85	0.96
3	1.86	0.99

Equations

$$\left\{\begin{array}[]{rcll}\partial_{t}u+\boldsymbol{\nabla}\cdot(\boldsymbol{\beta}u)-\boldsymbol{\nabla}\cdot(\mu\boldsymbol{\nabla}u)&=&g&\text{in }\Omega\times(0,T],\\ u(x,t)&=&\overline{u}(x,t)&\text{on }\partial\Omega\times(0,T],\\ u(x,0)&=&u_{0}(x)&x\in\Omega,\end{array}\right.$$

$$B_{i}^{0}(x)\doteq\left\{\begin{array}[]{ll}1&\text{ if }\ {\xi}_{i}\leq{\xi}<{\xi}_{i+1}\\ 0&\text{ otherwise}\end{array},\right.\qquad B_{i}^{k}(x)\doteq\dfrac{x-{\xi}_{i}}{{\xi}_{i+k}-{\xi}_{i}}B_{i}^{k-1}(x)+\dfrac{{\xi}_{i+k+1}-x}{{\xi}_{i+k+1}-{\xi}_{i+1}}B_{i+1}^{k-1}(x),$$

$$({\partial }_t{u}_h,v_h)+(\mathbf{\beta }\cdot \nabla u_h,v_h)+\mu (\nabla u_h,\nabla v_h)=(g,v_h),\forall v_h\in V_h,$$

$$\underset{\mathbf{j}\in {\mathcal{N}}_h^{\mathbf{i}}\backslash \{\mathbf{i}\}}{\min }{u}_{\mathbf{j}}\le u_{\mathbf{i}}\le \underset{\mathbf{j}\in {\mathcal{N}}_h^{\mathbf{i}}\backslash \{\mathbf{i}\}}{\max }{u}_{\mathbf{j}}.$$

$$\min \left(\underset{\mathbf{x}\in \partial \Omega ,t\in [0,T)}{\min }{\overline{u}}_h(\mathbf{x},t),\underset{\mathbf{x}\in \Omega }{\min }{u}_{h0}(\mathbf{x})\right)\le u_{\mathbf{i}}\le \max \left(\underset{\mathbf{x}\in \partial \Omega ,t\in [0,T)}{\max }{\overline{u}}_h(\mathbf{x},t),\underset{\mathbf{x}\in \Omega }{\max }{u}_{h0}(\mathbf{x})\right).$$

$$B_h(w_h;u_h,v_h)\doteq \sum _{\mathbf{i}\in {\mathcal{N}}_h}\sum _{\mathbf{j}\in {\mathcal{S}}_h^{\mathbf{i}}}{\nu }_{\mathbf{ij}}(w_h)v_{\mathbf{i}}{u}_{\mathbf{j}}\ell (\mathbf{i},\mathbf{j}),$$

$${\nu }_{\mathbf{ij}}(w_h)$$

$${\nu }_{\mathbf{i}\mathbf{i}}(w_h)$$

$$\left\llbracket\boldsymbol{\nabla}u_{h}\right\rrbracket_{\boldsymbol{ij}}\doteq\frac{u_{\boldsymbol{j}}-u_{\boldsymbol{i}}}{|\boldsymbol{r}_{\boldsymbol{ij}}|}+\frac{u_{\boldsymbol{j}}^{\rm sym}-u_{\boldsymbol{i}}}{|\boldsymbol{r}_{\boldsymbol{ij}}^{\rm sym}|},$$

$$\{\{|\nabla u_h\cdot {\hat{\mathbf{r}}}_{\mathbf{ij}}|\}{\}}_{\mathbf{ij}}\doteq \frac{1}{2}\left(\frac{|u_{\mathbf{j}}-u_{\mathbf{i}}|}{|{\mathbf{r}}_{\mathbf{ij}}|}+\frac{|u_{\mathbf{j}}^{\mathrm{sym}}-u_{\mathbf{i}}|}{|{\mathbf{r}}_{\mathbf{ij}}^{\mathrm{sym}}|}\right).$$

$$\alpha_{{\boldsymbol{i}}}(u_{h})\doteq\left\{\begin{array}[]{cc}\left[\dfrac{\left|\sum_{{\boldsymbol{j}}\in\mathcal{N}_{h}^{{\boldsymbol{i}}}}\left\llbracket\boldsymbol{\nabla}u_{h}\right\rrbracket_{\boldsymbol{ij}}\right|}{\sum_{{\boldsymbol{j}}\in\mathcal{N}_{h}^{{\boldsymbol{i}}}}2\{\kern-3.125pt\{\left|\boldsymbol{\nabla}u_{h}\cdot\hat{\boldsymbol{r}}_{\boldsymbol{ij}}\right|\}\kern-3.125pt\}_{\boldsymbol{ij}}}\right]^{q}&\text{if}\quad\sum_{{\boldsymbol{j}}\in\mathcal{N}_{h}^{{\boldsymbol{i}}}}\{\kern-3.125pt\{\left|\boldsymbol{\nabla}\cdot\hat{\boldsymbol{r}}_{\boldsymbol{ij}}\right|\}\kern-3.125pt\}_{\boldsymbol{ij}}\neq 0\\ 0&\text{otherwise}\end{array}\right..$$

$$({\partial }_t{u}_h,v_h)+(\mathbf{\beta }\cdot \nabla u_h,v_h)+\mu (\nabla u_h,\nabla v_h)+B_h(u_h;u_h,v_h)=(g,v_h),\forall v_h\in V_h,$$

$${\tilde{\mathsf{K}}}_{\mathbf{ij}}(u_h)u_{\mathbf{j}}={\mathsf{F}}_{\mathbf{i}},$$

$${\tilde{\mathsf{K}}}_{\mathbf{ij}}(u_h)\le 0,\forall \mathbf{j}\in {\mathcal{S}}_h^{\mathbf{i}}\backslash \{\mathbf{i}\},\sum _{\mathbf{j}\in {\mathcal{S}}_h^{\mathbf{i}}}{\tilde{\mathsf{K}}}_{\mathbf{ij}}(u_h)=0.$$

$$\sum _{\mathbf{j}\in {\mathcal{S}}_h^{\mathbf{i}}}{\tilde{\mathsf{K}}}_{\mathbf{ij}}(u_h)u_{\mathbf{j}}=0.$$

$$u_{\mathbf{i}}=\frac{{\sum }_{\mathbf{j}\in {\mathcal{S}}_h^{\mathbf{i}}\backslash \{\mathbf{i}\}}{\tilde{\mathsf{K}}}_{\mathbf{ij}}(u_h)u_{\mathbf{j}}}{{\tilde{\mathsf{K}}}_{\mathbf{i}\mathbf{i}}(u_h)}.$$

$$u_{\mathbf{i}}=\frac{{\tilde{\mathsf{K}}}_{\mathbf{i}\mathbf{k}}(u_h)u_{\mathbf{k}}}{{\tilde{\mathsf{K}}}_{\mathbf{i}\mathbf{i}}(u_h)}+\frac{{\sum }_{\mathbf{j}\in {\mathcal{S}}_h^{\mathbf{i}}\backslash \{\mathbf{i},\mathbf{k}\}}{\tilde{\mathsf{K}}}_{\mathbf{ij}}(u_h)u_{\mathbf{j}}}{{\tilde{\mathsf{K}}}_{\mathbf{i}\mathbf{i}}(u_h)}.$$

$$\left(1-\frac{{\tilde{\mathsf{K}}}_{\mathbf{i}\mathbf{k}}(u_h)}{{\tilde{\mathsf{K}}}_{\mathbf{i}\mathbf{i}}(u_h)}\right)u_{\mathbf{i}}=\frac{{\sum }_{\mathbf{j}\in {\mathcal{S}}_h^{\mathbf{i}}\backslash \{\mathbf{i},\mathbf{k}\}}{\tilde{\mathsf{K}}}_{\mathbf{ij}}(u_h)u_{\mathbf{j}}}{{\tilde{\mathsf{K}}}_{\mathbf{i}\mathbf{i}}(u_h)}.$$

$$u_{\mathbf{i}}=\frac{{\sum }_{\mathbf{j}\in {\mathcal{S}}_h^{\mathbf{i}}\backslash \{\mathbf{i},\mathbf{k}\}}{\tilde{\mathsf{K}}}_{\mathbf{ij}}(u_h)u_{\mathbf{j}}}{\left(1-\frac{{\tilde{\mathsf{K}}}_{\mathbf{i}\mathbf{k}}(u_h)}{{\tilde{\mathsf{K}}}_{\mathbf{i}\mathbf{i}}(u_h)}\right){\tilde{\mathsf{K}}}_{\mathbf{i}\mathbf{i}}(u_h)}=\frac{{\sum }_{\mathbf{j}\in {\mathcal{S}}_h^{\mathbf{i}}\backslash \{\mathbf{i},\mathbf{k}\}}{\tilde{\mathsf{K}}}_{\mathbf{ij}}(u_h)u_{\mathbf{j}}}{{\tilde{\mathsf{K}}}_{\mathbf{i}\mathbf{i}}(u_h)-{\tilde{\mathsf{K}}}_{\mathbf{i}\mathbf{k}}(u_h)}.$$

$$|w{|}_{\ell }\doteq \sqrt{\frac{1}{2}\sum _{\mathbf{i}\in {\mathcal{N}}_h}\sum _{j\in {\mathcal{S}}_h^{\mathbf{i}}}(w_{\mathbf{i}}-w_{\mathbf{j}}{)}^2}.$$

$$(\mathsf{B}(u)-\mathsf{B}(v),z)\le Cq(h^d+\Vert \mathbf{\beta }{\Vert }_{\infty }{h}^{d-1}\delta t+\mu h^{d-2}\delta t)|u-v{|}_{\ell }|z{|}_{\ell }$$

$${\mathsf{K}}_{\mathbf{ij}}=({\partial }_t{\phi }_{\mathbf{j}},{\phi }_{\mathbf{i}})+(\mathbf{\beta }\cdot \nabla {\phi }_{\mathbf{j}},{\phi }_{\mathbf{i}})+\mu (\nabla {\phi }_{\mathbf{j}},\nabla {\phi }_{\mathbf{i}}).$$

$${\mathsf{K}}_{\mathbf{ij}}$$

$$\le C_1{h}^d+\Vert \mathbf{\beta }{\Vert }_{\infty }{C}_2{h}^{d-1}\delta t+\mu C_3{h}^{d-2}\delta t,$$

$$({\partial }_t{u}_h^l,v_h)+(\mathbf{\beta }\cdot \nabla u_h^l,v_h)+\mu (\nabla u_h^l,\nabla v_h)+B_h(u_h^l;u_h^l,v_h)=(g,v_h),\forall v_h\in V_h^l,$$

$$\alpha_{\varepsilon_{h},i}(u_{h})\doteq\left[Z\left(\dfrac{\left|\sum_{{\boldsymbol{j}}\in\mathcal{N}_{h}^{{\boldsymbol{i}},s}}\left\llbracket\boldsymbol{\nabla}u_{h}\right\rrbracket_{\boldsymbol{ij}}\right|_{1,\varepsilon_{h}}+\gamma_{h}}{\sum_{{\boldsymbol{j}}\in\mathcal{N}_{h}^{{\boldsymbol{i}},s}}2\{\kern-3.125pt\{\left|\boldsymbol{\nabla}u_{h}\cdot\hat{\boldsymbol{r}}_{\boldsymbol{ij}}\right|_{2,\varepsilon_{h}}\}\kern-3.125pt\}_{\boldsymbol{ij}}+\gamma_{h}}\right)\right]^{q},$$

$${\left|x\right|}_{1,{\epsilon }_h}=\sqrt{x^2+{\epsilon }_h},{\left|x\right|}_{2,{\epsilon }_h}=\frac{x^2}{\sqrt{x^2+{\epsilon }_h}}.$$

$$Z\left(x\right)\doteq\left\{\begin{array}[]{ll}2x^{4}-5x^{3}+3x^{2}+x&x<1\\ 1&x\geq 1\end{array}\right.,$$

$$\max {}_{{\sigma }_h}\{x,y\}\doteq \frac{{\left|x-y\right|}_{1,{\sigma }_h}}{2}+\frac{x+y}{2}.$$

$$\tilde{\nu}_{\boldsymbol{ij}}(w_{h})\doteq\left\{\begin{array}[]{l}\max{}_{\sigma_{h}}\left\{\max{}_{\sigma_{h}}\left\{\alpha_{\varepsilon_{h},{\boldsymbol{i}}}(w_{h})\boldsymbol{\mathsf{K}}_{\boldsymbol{ij}},\alpha_{\varepsilon_{h},{\boldsymbol{j}}}(w_{h})\boldsymbol{\mathsf{K}}_{\boldsymbol{ji}}\right\},0\right\}\quad\text{for}\quad{\boldsymbol{j}}\neq{\boldsymbol{i}}\\ \displaystyle\sum_{{\boldsymbol{j}}\in\mathcal{S}_{h}^{{\boldsymbol{i}}}\backslash\{{\boldsymbol{i}}\}}\tilde{\nu}_{\boldsymbol{ij}}(w_{h})\end{array}\right.,$$

$${\tilde{B}}_h(w_h;u_h,v_h)\doteq \sum _{\mathbf{i}\in {\mathcal{N}}_h}\sum _{\mathbf{j}\in {\mathcal{S}}_h^{\mathbf{i}}}{\tilde{\nu }}_{ij}(w_h)v_{\mathbf{i}}{u}_{\mathbf{j}}\ell (\mathbf{i},\mathbf{j}).$$

Full text

† CIMNE – Centre Internacional de Mètodes Numèrics en Enginyeria, Esteve Terradas 5, 08860 Castelldefels, Spain ({jbonilla,sbadia}@cimne.upc.edu).

‡ Department of Civil and Environmental Engineering. Universitat Politècnica de Catalunya, Jordi Girona 1-3, Edifici C1, 08034 Barcelona, Spain.

Keywords: Isogeometric analysis, discrete maximum principle, monotonicity, high-order, space–time

1. 1 Introduction

2. 2 Preliminaries

3. 2.1 Convection-Diffusion problem

4. 2.2 Discretization

5. 2.3 Discrete problem

6. 2.4 Monotonicity properties

7. 3 Lipschitz-continuous nonlinear stabilization

8. 4 Time partitioned scheme

9. 5 Differentiable stabilization

10. 6 Numerical experiments

11. 6.1 1D Transient Diffusion

12. 6.2 Steady convection

13. 6.3 Nonlinear convergence

14. 6.4 1D Sharp layer propagation

15. 6.5 Boundary layer

16. 6.6 Three Body rotation

17. 7 Conclusions

FE finite element DOF degrees of freedom SSP strong stability preserving DMP discrete maximum principle MP maximum principle LED local extremum diminishing DLED discrete local extremum diminishing

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-Kutta 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 allows 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 $\Omega\times(0,T)\doteq\prod_{\alpha=1}^{d+1}(0,L_{\alpha})$ be a $(d+1)$-cube, where $d$ is the number of spatial dimensions. Then, the problem reads:

[TABLE]

where $\boldsymbol{\beta}$ is a divergence-free convection velocity, $\mu\geq 0$ is a scalar constant diffusion, and $g(\boldsymbol{x},t)$ is the body force. In the case of pure convection ($\mu=0$), boundary conditions are only imposed at the inflow $\Gamma_{\rm in}\doteq\{\boldsymbol{x}\in\partial\Omega:\boldsymbol{\beta}\cdot\boldsymbol{n}_{\partial\Omega}<0\}$, where $\boldsymbol{n}_{\partial\Omega}$ is a unit vector outward-pointing normal to the boundary. We also define the outflow boundary as $\Gamma_{\rm out}\doteq\partial\Omega\backslash\Gamma_{\rm in}$. Moreover, we will also consider the steady problem, which is obtained by dropping the time derivative term and the initial condition. 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 satisfies 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.

In order to avoid technicalities and facilitate the exposition of the stabilization method, we restrict this work to cubic domains. However, it is possible to work with complex geometries using standard procedures from isogeometric analysis [10]. E.g., a complex geometry would be divided in several parts, which would be mapped to multiple $d$-cubic patches. The stabilization method presented in this work is independent from this procedure.

2.2. Discretization

In this work, we consider a standard B-spline discretization with interpolative boundaries (see [10]). A spline of order $p$ in the variable $x$ is a piecewise polynomial function in $x$ of degree $p$. The values of $x$ in which different polynomials meet are called knots. Knots might be placed at the same location, i.e. can be repeated. When the knots are not repeated, the first $p-1$ derivatives of the spline are continuous. When a knot is repeated $r$ times, only the first $p-r$ derivatives are continuous across that knot. Knots are sorted in increasing order and collected in the so called knot vector $\{{\xi}_{1},{\xi}_{2},\ldots\}$. Given a knot vector, B-splines of order $p$ are basis functions for spline functions of the same order. B-splines are constructed in a recursive way using the Cox–de Boor formula:

[TABLE]

for $k=1,...,p$. By construction, $B_{i}^{p}(x)$ has compact support, is non-negative, and non-zero in $[{\xi}_{i},{\xi}_{i+p+1}]$. Notice that its support increases with the degree of the polynomial.

Let us consider the domain $[0,L]$ and the uniform partition into $m$ sub-intervals of size $h=L/m$. The open knot vector $\{{\xi}_{1},\ldots,{\xi}_{m+2p+1}\}$ is defined as follows. The first $p+1$ knots are located at zero, i.e., $\xi_{1}=\ldots=\xi_{p+1}=0$. The last $p+1$ knots are located at $L$, i.e., $\xi_{m+p+1}=\ldots=\xi_{m+2p+1}=L$. The interior points are equi-distributed, with $\xi_{i}=(i-p-1)h$, for $i=p+1,\ldots,m+p+1$. It leads to a basis $B_{i}^{p}(x)$ (for $i=1,\ldots,m+p$) for a space of splines in $[0,L]$ and a partition of unity, i.e., $\sum_{i=1}^{m}B_{i}^{p}(x)=1$ for $x\in[0,L]$. Any spline $v(x)$ of order $p$ in $[0,L]$ can uniquely be defined by the control points $(v_{1},\ldots,v_{m+p})\in\mathbb{R}^{m+d}$ as the linear combination of B-splines $v(x)=\sum_{i=1}^{m+p}B_{i}^{p}(x)v_{i}$. In one dimension, the basis functions obtained from an open knot vector are interpolatory at the extremes, i.e., $v(0)=v_{1}$ and $v(L)=v_{m+p+1}$ (see Fig. 1). For a first order polynomial $v$ in $[0,L]$, it holds $v(x)=\sum_{i=1}^{m}B_{i}^{p}(x)v(x_{i})$, where $x_{i}\doteq({\xi}_{i+1}+\ldots+{\xi}_{i+p})/p$ are called the Greville abscissae [11, 24].

Let us consider the number of partitions per dimension with $m_{\alpha}$, for $\alpha=1,\ldots,d+1$. We represent with $\mathcal{N}_{h}$ the set of multi-indices ${\boldsymbol{i}}\doteq(i_{1},...,i_{d+1})\in\mathbb{Z}^{d+1}$ with $i_{\alpha}\in\{1,\ldots,m_{\alpha}+p\}$. Every ${\boldsymbol{i}}\in\mathcal{N}_{h}$ can be expressed as $({{\boldsymbol{i}}_{\boldsymbol{x}}},{i_{t}})$, where ${{\boldsymbol{i}}_{\boldsymbol{x}}}$ is the spatial index and ${i_{t}}$ is the temporal index. The $(d+1)$-dimensional B-spline is defined as the tensor product of $d+1$ unidimensional B-splines ${B}_{\boldsymbol{i}}^{p}(\boldsymbol{x})\doteq B_{i_{1}}^{p}(x_{1})\times\cdots\times B_{i_{d+1}}^{p}(x_{d+1})$. Notice that a Greville abscissa in the case of a multidimensional spline reads $\boldsymbol{x}_{\boldsymbol{i}}=(x_{{\boldsymbol{i}}_{1}},...,x_{{\boldsymbol{i}}_{d+1}})$.

We define the space of splines $V_{h}\doteq\operatorname{span}\{{B}_{{\boldsymbol{i}}}^{p}(\boldsymbol{x})\ :\ {\boldsymbol{i}}\in\mathcal{N}_{h}\}$. We use the notation $\varphi_{{\boldsymbol{i}}}\equiv B_{\boldsymbol{i}}^{p}$. The order is omitted since it is assumed to be fixed. Thus, every spline $v_{h}\in V_{h}$ can be written as $v_{h}=\sum_{{\boldsymbol{i}}\in\mathcal{N}_{h}}\varphi_{{\boldsymbol{i}}}v_{\boldsymbol{i}}$. Furthermore, we define the following sets of indices, which are useful for the definition of the forthcoming schemes. The set of neighbors of ${\boldsymbol{i}}$ is defined as $\mathcal{N}_{h}^{{\boldsymbol{i}}}\doteq\{{\boldsymbol{j}}\in\mathcal{N}_{h}:|{\boldsymbol{i}}-{\boldsymbol{j}}|_{\infty}\leq 1\}$. We define as $\mathcal{S}_{h}^{{\boldsymbol{i}}}\doteq\{{\boldsymbol{j}}\in\mathcal{N}_{h}:|{\boldsymbol{i}}-{\boldsymbol{j}}|_{\infty}\leq p\}$ the set of indices whose associated shape functions intersect with the support of $\varphi_{{\boldsymbol{i}}}$.

We use standard notation for Sobolev spaces. The $L^{2}(\omega)$ scalar product is denoted by $(\cdot,\cdot)_{\omega}$ for $\omega\subset\Omega$. However, we omit the subscript for $\omega\equiv\Omega$. The $L^{2}(\Omega)$ norm is denoted by $\|\cdot\|$.

2.3. Discrete problem

The weak form of (1) using the Galerkin method reads: find $u_{h}\in V_{h}$ such that $u_{h}(\boldsymbol{x},t)=\overline{u}_{h}(\boldsymbol{x},t)$ on $\partial\Omega\times(0,T]$, $u_{h}(\boldsymbol{x},0)=u_{0h}(\boldsymbol{x})$ on $\Omega\times\{0\}$, and

[TABLE]

where $\overline{u}_{h}(t)$ and $u_{0h}$ are projections of $\overline{u}(t)$ and $u_{0}$ to $V_{h}$, respectively, such that the local DMP is satisfied (see Def. 2.2). Furthermore, we can rewrite the previous discrete problem in matrix form as $\boldsymbol{\mathsf{K}}_{\boldsymbol{ij}}u_{\boldsymbol{j}}=\boldsymbol{\mathsf{F}}_{\boldsymbol{i}}$, where $\boldsymbol{\mathsf{K}}_{\boldsymbol{ij}}\doteq(\partial_{t}\varphi_{{\boldsymbol{j}}},\varphi_{{\boldsymbol{i}}})+(\boldsymbol{\beta}\cdot\boldsymbol{\nabla}\varphi_{{\boldsymbol{j}}},\varphi_{{\boldsymbol{i}}})+\mu(\boldsymbol{\nabla}\varphi_{{\boldsymbol{j}}},\boldsymbol{\nabla}\varphi_{{\boldsymbol{i}}})$, and $\boldsymbol{\mathsf{F}}_{\boldsymbol{i}}\doteq(g,\varphi_{{\boldsymbol{i}}})$ for ${\boldsymbol{i}},{\boldsymbol{j}}\in\mathcal{N}_{h}$. Notice that we have not applied the boundary conditions yet. To apply boundary conditions the space of test functions is restricted to $v_{h}\in V_{h0}$, and the force vector is redefined as $\boldsymbol{\mathsf{F}}_{\boldsymbol{i}}\doteq(g,\varphi_{{\boldsymbol{i}}})-(\partial_{t}\overline{u}_{h},\varphi_{{\boldsymbol{i}}})-{(\boldsymbol{\beta}\cdot\boldsymbol{\nabla}\overline{u}_{h},\varphi_{{\boldsymbol{i}}})}-\mu(\boldsymbol{\nabla}\overline{u}_{h},\boldsymbol{\nabla}\varphi_{{\boldsymbol{i}}})$.

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 $u_{h}$ 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 $u_{h}\in V_{h}$ has a local discrete minimum (resp. maximum) on ${\boldsymbol{i}}\in\mathcal{N}_{h}$, if $u_{\boldsymbol{i}}\leq u_{\boldsymbol{j}}$ (resp. $u_{\boldsymbol{i}}\geq u_{\boldsymbol{j}}$) $\forall{\boldsymbol{j}}\in\mathcal{N}_{h}^{{\boldsymbol{i}}}$.

For problems that satisfy a maximum principle, e.g., problem (3) with $g=0$, 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 $u_{h}\in V_{h}$ satisfies the local discrete maximum principle if for every ${\boldsymbol{i}}\in\mathcal{N}_{h}$

[TABLE]

Definition 2.3** (Global space–time DMP).**

A solution $u_{h}\in V_{h}$ satisfies the global discrete maximum principle if the global extrema are located at boundary conditions, i.e., for every ${\boldsymbol{i}}\in\mathcal{N}_{h}$

[TABLE]

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, $\boldsymbol{\mathsf{B}}_{\boldsymbol{ij}}(u_{h})$, is said to be linearity-preserving if, for a solution that is linear in all directions in the neighborhood of $\boldsymbol{x}_{\boldsymbol{i}}$, then the stabilization term becomes null, i.e., $\boldsymbol{\mathsf{B}}_{\boldsymbol{ij}}(u_{h})=0$ if $u_{h}(\boldsymbol{x})\in\mathcal{P}_{1}(\Omega_{\boldsymbol{i}})$ where $\Omega_{\boldsymbol{i}}$ is the convex hull defined by the set of neighboring Greville abscissae $\{\boldsymbol{x}_{\boldsymbol{j}}\}_{{\boldsymbol{j}}\in\mathcal{N}_{h}^{\boldsymbol{k}},\ \boldsymbol{k}\in\mathcal{S}_{h}^{{\boldsymbol{i}}}}$.

3. Lipschitz-continuous nonlinear stabilization

In this section we define a nonlinear stabilization operator, $B_{h}(w_{h};u_{h},v_{h})$, to be added to the discrete problem (3), such that it satisfies at least the global DMP in Def. 2.3. Let us define $\boldsymbol{\mathsf{B}}_{\boldsymbol{ij}}(u_{h})\doteq B_{h}(u_{h};\varphi_{{\boldsymbol{j}}},\varphi_{{\boldsymbol{i}}})$. We also enforce that, for any $u_{h}\in V_{h}$, $\boldsymbol{\mathsf{B}}_{\boldsymbol{ij}}(u_{h})$

1. (1)

has compact support: $\boldsymbol{\mathsf{B}}_{\boldsymbol{ij}}(u_{h})=0$ if ${\boldsymbol{j}}\not\in\mathcal{S}_{h}^{{\boldsymbol{i}}}$, 2. (2)

is symmetric: $\boldsymbol{\mathsf{B}}_{\boldsymbol{ij}}(u_{h})=\boldsymbol{\mathsf{B}}_{\boldsymbol{ji}}(u_{h})$, 3. (3)

is conservative: $\sum_{{\boldsymbol{j}}\in\mathcal{S}_{h}^{{\boldsymbol{i}}}\backslash\{{\boldsymbol{i}}\}}\boldsymbol{\mathsf{B}}_{\boldsymbol{ij}}(u_{h})=-\boldsymbol{\mathsf{B}}_{{\boldsymbol{i}}{\boldsymbol{i}}}(u_{h})$.

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

[TABLE]

for any $w_{h},\ u_{h},\ v_{h}\in V_{h}$. Here, $\ell$ is the graph-Laplacian operator defined as $\ell({\boldsymbol{i}},{\boldsymbol{j}})=2\delta_{\boldsymbol{ij}}-1$ (see [13, 2]), and $\nu_{\boldsymbol{ij}}(w_{h})$ is the artificial diffusion defined as

[TABLE]

We denote by $\alpha(w_{h})$ 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 $\nu_{\boldsymbol{ij}}$. A shock detector must be a positive real number, which takes value 1 when $u_{h}(\boldsymbol{x}_{\boldsymbol{i}})$ is an inadmissible value of $u_{h}$ (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 $u_{h}\in\mathcal{P}_{1}(\Omega\times(0,T])$. In this section, we propose an isotropic approach for $\alpha_{{\boldsymbol{i}}}(w_{h})$, which consists in using the shock detector in [2] in all directions (including time).

In order to introduce the shock detector, let us recall some useful notation from [2]. Let $\boldsymbol{r}_{{\boldsymbol{ij}}}=\boldsymbol{x}_{\boldsymbol{j}}-\boldsymbol{x}_{\boldsymbol{i}}$ be the vector pointing from Greville abscissae $\boldsymbol{x}_{\boldsymbol{i}}$ to $\boldsymbol{x}_{\boldsymbol{j}}$ with ${\boldsymbol{i}},{\boldsymbol{j}}\in\mathcal{N}_{h}$ and $\hat{\boldsymbol{r}}_{{\boldsymbol{ij}}}\doteq\frac{\boldsymbol{r}_{{\boldsymbol{ij}}}}{|\boldsymbol{r}_{{\boldsymbol{ij}}}|}$. Let us take the set of Greville abscissae $\boldsymbol{x}_{\boldsymbol{j}}$ for ${\boldsymbol{j}}\in\mathcal{N}_{h}^{i}\backslash\{{\boldsymbol{i}}\}$ as vertices of a polytope in $d+1$ dimensions. In particular, let us name this polytope $\mathcal{Q}_{\boldsymbol{i}}$. Let $\boldsymbol{x}_{\boldsymbol{ij}}^{\rm sym}$ be the point at the intersection between $\partial\mathcal{Q}_{\boldsymbol{i}}$ and the line that passes through $\boldsymbol{x}_{\boldsymbol{i}}$ and $\boldsymbol{x}_{\boldsymbol{j}}$ that is not $\boldsymbol{x}_{\boldsymbol{j}}$ (see Fig. 2). The set of all $\boldsymbol{x}_{\boldsymbol{ij}}^{\rm sym}$ for all ${\boldsymbol{j}}\in\mathcal{N}_{h}^{{\boldsymbol{i}}}\backslash\{{\boldsymbol{i}}\}$ is represented with $\mathcal{N}^{{{\boldsymbol{i}}},{\rm sym}}_{h}$. We define $\boldsymbol{r}_{\boldsymbol{ij}}^{\rm sym}\doteq\boldsymbol{x}_{\boldsymbol{ij}}^{\rm sym}-\boldsymbol{x}_{\boldsymbol{i}}$. Given $\boldsymbol{x}^{\rm sym}_{\boldsymbol{ij}}$ in two dimensions, let as call $\boldsymbol{a}$ and $\boldsymbol{b}$ the indices of the vertices such that they define the edge in $\partial\mathcal{Q}_{\boldsymbol{i}}$ that contains $\boldsymbol{x}_{\boldsymbol{ij}}^{\rm sym}$. We define $u_{\boldsymbol{j}}^{\rm sym}$ as the linear interpolation of $u_{\boldsymbol{a}}$ and $u_{\boldsymbol{b}}$ at $\boldsymbol{x}_{\boldsymbol{ij}}$, i.e. $u_{\boldsymbol{ij}}^{\rm sym}\doteq u_{\boldsymbol{a}}\frac{\boldsymbol{x}_{\boldsymbol{b}}-\boldsymbol{x}_{\boldsymbol{ij}}^{\rm sym}}{\boldsymbol{x}_{\boldsymbol{b}}-\boldsymbol{x}_{\boldsymbol{a}}}+u_{\boldsymbol{b}}\frac{\boldsymbol{x}_{\boldsymbol{ij}}^{\rm sym}-\boldsymbol{x}_{\boldsymbol{a}}}{\boldsymbol{x}_{\boldsymbol{b}}-\boldsymbol{x}_{\boldsymbol{a}}}$. For higher dimensions, $u_{\boldsymbol{ij}}^{\rm sym}$ is defined analogously. Given the facet of $\partial\mathcal{Q}_{i}$ where $\boldsymbol{x}_{\boldsymbol{ij}}^{\rm sym}$ lies, $u_{\boldsymbol{ij}}^{\rm sym}$ is the linear interpolation at $\boldsymbol{x}_{\boldsymbol{ij}}^{\rm sym}$ 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 $u_{h}\in\mathcal{P}_{1}$, $u_{h}(\boldsymbol{x}_{\boldsymbol{i}})=u_{\boldsymbol{i}}$. Therefore, one can construct easily linear approximations of the unknown gradients that are exact for $u_{h}\in\mathcal{P}_{1}$. Furthermore, one can define the jump and the mean of a linear approximation of the unknown gradient at Greville abscissa $\boldsymbol{x}_{\boldsymbol{i}}$ in direction $\boldsymbol{r}_{\boldsymbol{ij}}$ as

[TABLE]

[TABLE]

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

[TABLE]

From Lm. 3.1 in [2] we know that (11) is valued between 0 and 1, and it is only equal to one if $u_{h}(\boldsymbol{x}_{i})$ 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 $u_{h}\in\mathcal{P}_{1}$, 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 $u_{h}\in V_{h}$ such that $u_{h}=\overline{u}_{h}$ on $\partial\Omega$, $u_{h}=u_{0h}$ at $t=0$, and

[TABLE]

which in turn can be expressed in matrix form as

[TABLE]

where $\tilde{\boldsymbol{\mathsf{K}}}_{\boldsymbol{ij}}(u_{h})\doteq\boldsymbol{\mathsf{K}}_{\boldsymbol{ij}}+\boldsymbol{\mathsf{B}}_{\boldsymbol{ij}}(u_{h})$ for ${\boldsymbol{i}},{\boldsymbol{j}}\in\mathcal{N}_{h}$.

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 $g=0$ and, for every control point ${\boldsymbol{i}}\in\mathcal{N}_{h}$ such that $u_{\boldsymbol{i}}$ is a local discrete extremum, it holds:

[TABLE]

Moreover, the resulting scheme is linearity-preserving as defined in Def. 2.4, i.e. $\boldsymbol{\mathsf{B}}_{ij}(u_{h})=0$ for $u_{h}\in\mathcal{P}_{1}$.

Proof.

Let us assume that $u_{\boldsymbol{i}}$ is a discrete maximum. Then, (13), for $g=0$ and before applying boundary conditions, reads

[TABLE]

Therefore, $u_{\boldsymbol{i}}$ can be computed as

[TABLE]

Since $u_{\boldsymbol{i}}$ is an extremum, which implies $\alpha_{{\boldsymbol{i}}}=1$, the stabilization term ensures (14) by construction. Hence, the coefficients that multiply $u_{\boldsymbol{j}}$ are in $[0,1]$, and the sum of all these coefficients add up to one. Therefore, $u_{\boldsymbol{i}}$ is a convex combination of its neighbors (including boundary conditions $\overline{u}$). Since $u_{\boldsymbol{i}}$ is a maximum and a convex combination of its neighbors, then $u_{\boldsymbol{k}}=u_{\boldsymbol{i}}$ for some $\boldsymbol{k}\in\mathcal{S}_{h}^{i}$. In that case, we can write

[TABLE]

Since $u_{\boldsymbol{k}}=u_{\boldsymbol{i}}$,

[TABLE]

Therefore, it can also be proved that $u_{\boldsymbol{i}}$ is a convex combination of all its neighbors but $u_{\boldsymbol{k}}$,

[TABLE]

Proceeding analogously, one can also prove that $u_{\boldsymbol{k}}$ is a convex combination of all its neighbors but $u_{\boldsymbol{i}}$. Hence, we know that the value of $u_{\boldsymbol{i}}=u_{\boldsymbol{k}}$ is bounded by all their neighbors. At this point, the same reasoning can be applied to any of their neighbors. Thus, 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], $\alpha_{{\boldsymbol{j}}}=0$ for any ${\boldsymbol{j}}\in\mathcal{S}_{h}^{{\boldsymbol{i}}}$ if $u_{h}\in\mathcal{P}_{1}(\Omega_{\boldsymbol{i}})$ where $\Omega_{\boldsymbol{i}}$ is the convex hull defined by the set of neighboring Greville abscissae $\{\boldsymbol{x}_{\boldsymbol{j}}\}_{{\boldsymbol{j}}\in\mathcal{N}_{h}^{\boldsymbol{k}},\ \boldsymbol{k}\in\mathcal{S}_{h}^{{\boldsymbol{i}}}}$. By definition, the stabilization term also vanishes if $\alpha_{{\boldsymbol{j}}}=0$ for ${\boldsymbol{j}}\in\mathcal{S}_{h}^{{\boldsymbol{i}}}$ (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 $q\rightarrow\infty$.

Proof.

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

Finally, we prove Lipschitz continuity of the stabilization term. In this particular case, the proof follows the same reasoning as in [2]. Let us recall the definition of the semi-norm generated by the graph-Laplacian required to show Lipschitz continuity,

[TABLE]

In addition, we will also need the $L^{2}(\Omega)$ norm denoted as $\|\cdot\|$ and the $L^{\infty}(\Omega)$ expressed as $\|\cdot\|_{\infty}$. We do not include all details for the sake of conciseness and refer the reader to the previously cited work.

Theorem 3.3**.**

Let $V_{h}^{\rm adm}\subset V_{h}$ be the subspace of functions that satisfy the global DMP in Def. 2.3, then $\boldsymbol{\mathsf{B}}(\cdot)$ with the shock detector (11) is Lipschitz continuous in $V_{h}^{\rm adm}$ for $u_{h}\in V_{h}$ and bounded $q$, since

[TABLE]

is satisfied.

Proof.

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

[TABLE]

In this case, using Cauchy-Schwarz inequality, the inverse inequality $\|\boldsymbol{\nabla}v_{h}\|\leq Ch^{-1}\|v_{h}\|$ for $v_{h}\in V_{h}$, and $\|\varphi_{i}\|\leq Ch^{d/2}$, we get

[TABLE]

where $d$ is the number of spatial dimensions, $h$ is the distance between knots for the spatial directions and $\delta t$ is the distance between knots for the time direction. ∎

Notice that whereas $C$ is uniform with respect to the mesh size, it can depend on the polynomial order of the discretization.

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 $\eqref{continuous-problem}$ set in $\Omega\times(0,T]$, will be decomposed in $\Omega\times(t_{l},t_{l+1}]$ for $0\leq l\leq n_{t}-1$, with $t_{0}=0$ and $t_{n_{t}}=T$. We define the length of the subdomain as $\Delta t\doteq t_{l+1}-t_{l}$, and restrict its possible values to $\Delta t=n\,p\,\delta t$ for some $n\in\mathbb{N}$, where $p$ is the order of the spline space, and $\delta t$ 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 $\Delta t$, $t_{l}$ 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 $\Omega\times(0,T]$, we discretize it as described in Sect. 2.2, resulting in a spline space $V_{h}$. Then in order to reduce the coupling between partitions, we insert $p$ knots at $t_{l}$. The resulting spaces at each subdomain, say $V_{h}^{l}$, 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 $l$ will affect the solution at $l+1$, but not the opposite. Therefore, we impose that the initial conditions at subdomain $l+1$ are equal to the solution at the final time of subdomain $l$, i.e. $u_{h}^{l+1}(t_{l})\doteq u_{h}^{l}(t_{l})$. After imposing this restriction, the complete approximation space, $\tilde{V_{h}}$, is $C^{0}$, 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 $l=1,...,n_{t}$; find $u_{h}^{l}\in V_{h}^{l}$ such that $u_{h}^{l}=\overline{u}_{h}$ on $\partial\Omega$, $u_{h}^{l}(\boldsymbol{x},t_{l})=u_{h}^{l-1}(\boldsymbol{x},t_{l})$ with $u_{h}^{0}(\boldsymbol{x},t_{0})=u_{0h}$, and

[TABLE]

Due to the partition $u_{h}$ will only be piecewise continuous in time. Let us prove now that the scheme still satisfies the global DMP.

Lemma 4.1**.**

The solution of problem (25), using the shock detector defined in (11), satisfies the global DMP (Def. 2.3) if $g=0$ in $\Omega\times(0,T]$ and, for every control point ${\boldsymbol{i}}\in\mathcal{N}_{h}$ such that $u_{\boldsymbol{i}}$ is a local discrete extremum, conditions (14) hold.

Proof.

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 $t_{1}$, $u_{h}^{1}(\boldsymbol{x},t_{1})$, is bounded by the initial and boundary conditions. Since $u_{h}^{1}(\boldsymbol{x},t_{1})$ 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 $u_{h}^{1}(\boldsymbol{x},t_{1})$, 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

[TABLE]

where $\gamma_{h}>0$ is a parameter to prevent division by zero, and the regularized absolute values by

[TABLE]

Notice that $\left|x\right|_{2,\varepsilon_{h}}\leq|x|\leq\left|x\right|_{1,\varepsilon_{h}}$. 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 $Z\left(x\right)$, which reads

[TABLE]

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 $[0,1]$ and it is equal to 1 if $u_{\boldsymbol{i}}$ 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 $\max{}_{\sigma_{h}}\{\cdot,\cdot\}$ as

[TABLE]

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

[TABLE]

and the stabilization operator reads

[TABLE]

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

[TABLE]

where $d$ is the spatial dimension of the problem, $L$ is a characteristic length, and $\sigma,\ \varepsilon,$ and $\gamma$ 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 $t\in(0,1]$ and $x\in\Omega\doteq(0,1)$,

[TABLE]

where $f\doteq 2(6x^{2}-6x+1)(t(t-1))^{2}+2t(t-1)(2t-1)(x(x-1))^{2}$. This problem has $u=(x(x-1))^{2}\ (t(t-1))^{2}$ 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 ${\delta t=\{0.2,0.1,0.05,0.025,0.0125\}}$ for the first order discretization. In spatial directions, the distance is small enough $(h=1/400)$ to prevent that spatial discretization errors affect the analysis. Second and third order discretizations are obtained using the following $k$-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 $k$-refinements to $p=2$ and $p=3$ discretizations, for an interior subset of the discretization. Henceforth, we will use this kind of $k$-refinement in order to increase the discretization order.

We measure the relative $L^{2}$ norm and $H^{1}$ semi-norm of error in the whole space–time domain, and compute the resulting convergence rate. Errors in $L^{2}$ norm and $H^{1}$ 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.

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 $x\in\Omega\doteq[0,1]^{2}$,

[TABLE]

where $u_{D}=\sin\left(2\pi\left(x-\frac{y}{\tan\theta}\right)\right)$, $\boldsymbol{\beta}=(\cos\theta,\sin\theta)$, and $\theta=\pi/3$. The analytical solution of the above problem reads $u=\sin\left(2\pi\left(x-\frac{y}{\tan\pi/3}\right)\right)$. The convergence analysis is performed for first, second, and third order discretizations, i.e. $p=\{1,2,3\}$. We use a standard nonlinear solver (see [3] for details), with a nonlinear tolerance $\frac{u^{k+1}-u^{k}}{u^{k}}<10^{-6}$. The following stabilization parameters have been used: $q=10$, $\varepsilon=10^{-5}$, $\sigma=10^{-6}$, $\gamma=10^{-10}$. 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 $L^{2}$ error norm and first in the $H^{1}$ error semi-norm. It is known that the stabilized scheme should recover second order convergence for $p=1$. 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 $k$-refinement previously defined.

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.

[TABLE]

where $\boldsymbol{\beta}\doteq 1$, $\Omega\doteq(0,1]$, $T=0.5$, and $u_{0}\doteq 1-H_{0}(x-0.25)$, where $H_{0}$ is the well-known zero-centered Heaviside function. First and second discretization orders are used in a coarse mesh of $25\times 25$ control points. To obtain the second order mesh, we perform the $k$-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: $\gamma=10^{-10}$, $\sigma=\zeta$, $\varepsilon=\zeta^{2}$, where $\zeta=\{10^{-1},10^{-2},10^{-3},10^{-4}\}$. Furthermore, the effect is also compared as $q$ is incremented, particularly for $q=\{1,2,5,10\}$. 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, and the nonlinear tolerance is set to $10^{-5}$.

Fig. 6 and 7 show the results for first and second order discretizations, respectively. In general terms, as $q$ is increased or $\zeta$ 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.

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 (35). 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 $10^{-5}$. The regularization parameters used are $q=10$, $\varepsilon=10^{-8}$, $\sigma=10^{-6}$, and $\gamma=10^{-10}$. 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 $k$-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. 8(a) the solution at $t=0.5$ is shown for different orders and fixed number of control points, and using the $k$-refinement in Fig. 8(b). We observe that for non-smooth solutions, fixing the number of control points and increasing the order does not improve the results. This is a consequence of the underlying discretization properties. The support of the shape functions becomes larger as the order is increased. Therefore, nonsmooth solutions become slightly more smeared. In the case of Fig. 8(b), as expected, we observe better approximations as the order is increased using the $k$-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. 9(a) as we increase the discretization order. In a similar trend, we observe less improvement in Fig. 9(b) than in Fig. 8(b). We attribute this degradation to the time partitions, which becomes more evident as the subdomains are smaller. In particular, at the boundary of each partition the method might slightly increase the amount of diffusion introduced. At these boundaries, the shock detector relies on a smaller domain to determine if the DMP is satisfied. Therefore, it is more likely to introduce more diffusion. On the other hand, the partition itself modifies the scheme introducing some error as shown in Sect. 6.1.

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 $\Omega\doteq[0,1]^{2}$

[TABLE]

where $\boldsymbol{\beta}=(\cos\theta,\sin\theta)$, $\theta=-\pi/3$, and the boundary conditions are defined as

[TABLE]

For this test, we use the following settings: a nonlinear tolerance of $10^{-8}$, $q=2$, $\varepsilon=10^{-8}$, $\sigma=10^{-6}$, and $\gamma=10^{-10}$. In Fig. 10(a), the solution for $p=4$ is depicted. The converged solution does not exhibit any oscillation. Very sharp layers are obtained for this parameter setting. In Fig. 10(b), we show the profile of the solution at $y=0.1$ for different orders. In this case, we start with a discretization of 50 control points per direction. Then, we increase the order using the $k$-refinement used previously. As previously observed for transient problems, we observe an improvement of the solution as the order is increased.

6.6. Three Body rotation

Finally, we solve the transient pure convection problem (35) in $\Omega\times(0,1]$ for $\Omega=[0,1]^{2}$, with $\boldsymbol{\beta}=(-2\pi(y-0.5),2\pi(x-0.5))$. Initial conditions are given in [17]. Its interpolation in a first order $200\times 200$ control point mesh is depicted in Fig. 11(a). The analytical solution of this problem is simply the translation of the profiles in the direction of the convection. In particular, for $t=1$, 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 $t=1$ and comparing the results with the initial conditions.

The solution is computed using scheme (12) in combination with the shock detector in (26). We use the following parameters for the stabilization: $q=10$, $\sigma=10^{-6}$, $\varepsilon=10^{-8}$, and $\gamma=10^{-10}$. Different meshes, time partitions, and discretization orders are used in this experiment. We start with a linear discretization of $100\times 100$ control points in space, and 500 in time divided in 125 subdomains. Then, we increase the discretization order to $p=2$ using the $k$-refinement. In order to compare first and second order discretizations, but using a similar number of control points we use a discretization with $200\times 200$ 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 $100\times 100\times 500$ 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.

Fig. 13 shows the solutions for $100\times 100$ meshes, and 125 subdomains in time, whereas Fig. 14 shows 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 13(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 $p=2$. 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. 14(a) and 13(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 $s\doteq\{(x,y):\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}=0.25\}$ is depicted.

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 $hp$-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.

Acknowledgments

J. Bonilla gratefully acknowledges the support received from ”la Caixa” Foundation through its PhD scholarship program (LCF/BQ/DE15/10360010). 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.