{'fraction_non_alphanumeric': 0.04275680007542545, 'fraction_numerical': 0.012987319096780275, 'mean_word_length': 4.21664822328784, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Monotonicity-preserving finite element schemes with adaptive mesh refinement for hyperbolic problems

This work is focused on the extension and assessment of the monotonicity-preserving scheme in [3] and the local bounds preserving scheme in [5] to hierarchical octree adaptive mesh refinement (AMR). Whereas the former can readily be used on this kind of meshes, the latter requires some modifications. A key question that we want to answer in this work is whether to move from a linear to a nonlinear stabilization mechanism pays the price when combined with shock-adapted meshes. Whereas nonlinear (or shock-capturing) stabilization leads to improved accuracy compared to linear schemes, it also negatively hinders nonlinear convergence, increasing computational cost. We compare linear and nonlinear schemes in terms of the required computational time versus accuracy for several steady benchmark problems. Numerical results indicate that, in general, nonlinear schemes can be cost-effective for sufficiently refined meshes. Besides, it is also observed that it is better to refine further around shocks rather than use sharper shock capturing terms, which usually yield stiffer nonlinear problems. In addition, a new refinement criterion has been proposed. The proposed criterion is based on the graph Laplacian used in the definition of the stabilization method. Numerical results show that this shock detector performs better than the well-known Kelly estimator for problems with shocks or discontinuities.


Introduction
Natural phenomena can develop shock waves in different scenarios. A classical example is the shock wave generated by an object traveling faster than sound. The numerical modeling of problems with shocks is still a challenge, especially when the admissible physical solution has some physical constraints, e.g., positivity or non-negativity, that must be preserved at the discrete level to have well-posedness; E.g., the fluid density and temperature are positive quantities in a compressible flow.
Several numerical schemes have been proposed so far to approximate this kind of problems by combining finite volume methods (FVM) or discontinuous Galerkin (dG) finite elements (FEs) for space discretization with explicit time integrators (see [19,23,40,56]). Explicit time integrators are only stable under a Courant-Friedrichs-Levy (CFL) restriction over the time step size, which implies to capture all time scales. Thus, explicit methods are not suitable for problems in which the smallest time scales are not of interest. For instance, the fastest time scales at a confined plasma in a nuclear fusion reactor are not of engineering interest whereas explicit time integration is unaffordable in practical simulations [31].
Implicit monotonicity-preserving (or at least positivity-preserving) methods are still scarce. As proved by Godunov [24], linear monotonicity-preserving schemes can be at most first-order accurate. For scalar problems (and under some mesh restrictions), Burman and Ern [18], Barrenechea and co-workers [12,13], Kuzmin and co-workers [33,36,43], and Badia and Hierro [6,7] have proposed nonlinear schemes that preserve monotonicity and can presumably attain higher order accuracy. 1 However, these properties comes at the cost of solving a very stiff nonlinear problem [32]. The authors [3,4] have proposed differentiable schemes that improve the nonlinear convergence behavior of previous methods.
For hyperbolic systems of equations, numerical methods are even less well developed. For explicit time integration, Guermond and Popov [27] have recently proposed a continuous Galerkin (cG) FE scheme that preserves positivity of density and energy under certain CFL-like condition. Unfortunately, these ideas cannot be easily extended to implicit time integration and we are not aware of any implicit method that theoretically satisfies such properties. Kuzmin and co-workers [36,41,46,47] have proposed various schemes based on flux corrected transport (FCT) [44] that are experimentally robust, but lack of a theoretical analysis. Besides, this strategy also yields very stiff nonlinear problems. Differentiable schemes for compressible flows have been proposed in [5] to alleviate (but not eliminate) this problem.
Shocks are non-smooth and localized and thus suitable for AMR [21,57]. AMR allows one to increase the mesh resolution only in the vicinity of shocks or discontinuities. In brief, the AMR process can be divided into two main ingredients. On the one hand, to estimate the error at each element. On the other hand, to decide which elements need to be refined or coarsened. This iterative process provides a mesh locally adapted to the features of the problem at hand. As a result, it is a nonlinear approximation scheme which tries to minimize the error for a target computational cost. If performed optimally, AMR exhibits exponential convergence even for solutions with limited regularity [21].
In this context, a key question is whether it is computationally more effective to consider a nonlinear high-order scheme (with the nonlinear convergence issues) or a cheaper linear (first-order) scheme in a more refined mesh. The motivation of this work is to shed light on this issue. First, we adapt the schemes developed in [3,5] to hierarchical octree AMR [9,54]. Next, we propose a refinement criterium that relies on information already present in the stabilization technique; nonlinear stabilization methods include a shock detector to activate the artificial diffusion only close to discontinuities. We propose to use a modification of the shock detector in [3] to drive the AMR process.
This paper is structured as follows. First, we introduce the problem, its discretization, and monotonicity properties for scalar problems and hyperbolic systems in Sect. 2. Then, the stabilization 1 In this work, schemes with nonlinear stabilization are also referred to as high-order and linear stabilization schemes as low or first-order.
where u = {u β } m β=1 are m ≥ 1 conserved variables, f is the physical flux,ū β (x, t) are the boundary values for the βth-component of u, u 0 (x) are the initial conditions, and g(x, t) is a function defining the body forces. Note that the flux, f : We denote by f : R m → R m×m×d the flux Jacobian. Let n ∈ R d be any direction vector. Since the system is hyperbolic, the flux Jacobian in any direction is diagonalizable and has only real eigenvalues, i.e., These eigenvalues might have different multiplicities and different signs. Hence, for a given direction n, each characteristic variable might be convected forward (along n) or backwards (along −n). Therefore, it is convenient to define inflow and outflow boundaries for each component.
where n ∂Ω is the unit outward normal to the boundary and λ β is the βth-eigenvalue of the flux Jacobian. We define the outflow boundary as Γ β out . = ∂Ω\Γ β in . We refer the reader to [23,28,56] for a detailed discussion on boundary conditions for hyperbolic problems. In the present study, we will also consider the steady counterpart of (1), which is obtained by dropping the time derivative term and the initial conditions.
In this work, we work with both scalar convection equations and Euler equations. Taking m = 1 and f (u) . = vu with v a divergence-free convection field, we recover the well known scalar transport problem. On the other hand, Euler equations for ideal gases are recovered by defining m = d + 2 and where ρ is the density, E is the total energy, p is the pressure, m = {m 1 , . . . , m d }, where m i = ρv i , is the momentum, v = {v 1 , . . . , v d } is the velocity, b = {b 1 , . . . , b d } are the body forces, r is an energy source term per unit mass, and I is the identity matrix of dimension d × d. In addition, the system is equipped with the ideal gas equation of state p = (γ − 1)ρı, where ı = E − 1 2 v 2 is the internal energy and γ is the adiabatic index.

2.2.
Discretization. The discretization used in this work is able to adapt its local size to the features of the problem at hand. In particular, it is a hierarchically refined octree-based hexahedral mesh [54]. This type of discretizations are constructed hierarchically. At every step of the refinement process, marked cells are refined into four (eight) cells in 2D (3D). The adaptation of the mesh to the problem at hand is achieved by only marking for refining a targeted amount of cells. This results in a mesh with different refinement levels at different regions. Hanging nodes appear at the interface between cells at different refinement levels. These are nodes that only belong to the cells at a higher refinement level (see Fig. 1). In our case, the meshes used are 2:1 balanced. This restriction implies that there can only be a difference of one refinement level between neighboring cells. This restriction is a trade-off between implementation complexity and performance gain that has been adopted by many AMR codes [54].
Hanging nodes need to be treated carefully in the case of working with conforming FE discretizations. Otherwise, associating a regular degree of freedom (DOF) to a hanging node may lead to discontinuities in the approximated solution. To preserve continuity of the FE space, hanging DOFs values are not included in the assembled system of equations but obtained by interpolating the values of the neighboring regular DOFs. For more details in the definitions of these constraints we refer the reader to [8][9][10].
Let T h be a hierarchical octree-based partition of Ω. Consider a Lagrangian (nodal) FE space on top of this mesh. The set of all nodes in the FE space is represented withÑ h . For every node i ∈Ñ h , x i stands for the node coordinates. We can splitÑ h into two subsets, namely the set of hanging nodes N hg h and the set of conforming nodes N h = card(N h ) the total number of conforming nodes. The set of nodes belonging to a particular element K ∈ T h is defined as Moreover, Ω i is the macroelement composed by the union of elements that contain node i, i.e., Ω i . = K∈T h , xi∈K K. To simplify the discussion below, we abuse notation and use i for both the node and its associated index.
We restrict the present work to first order FEs and define the FE space as follows. We define is the space of polynomials of partial degree less than or equal to one. Furthermore, we define the space V h0 .
The functions v h ∈ V h can be constructed as a linear combination of the basis {ϕ i } i∈Ñ h and nodal values v i , where ϕ i is the shape function associated to the node i. Hence, v h = i∈Ñ h ϕ i v i .
We use standard notation for Sobolev spaces. The L 2 (ω) scalar product is denoted by (·, ·) ω for ω ⊂ Ω. However, we omit the subscript for ω ≡ Ω. The L 2 norm is denoted by · .
The method of lines is applied in combination with the FE spaces described above for the spatial discretization. We approximate the solution u ≈ u h = i∈Ñ h ϕ i u i , and the fluxes as f ≈ f h = i∈Ñ h ϕ i f (u i ). For simplicity in the exposition, we use the Backward Euler (BE) scheme for the time discretization; higher order time discretizations can be achieved using strong stability preserving (SSP)-Runge Kutta (RK) methods (see [25]). In the latter case, a CFL-like condition must be satisfied to enjoy the monotonicity properties in Sect. 2.3 (see [35,39]).
The semi-discrete Galerkin FE approximation of the weak form of (1) reads: where u β h and u 0h are admissible FE approximations of the boundary and initial conditions u β and u 0 . In this context, we consider admissible any approximation that satisfies the maximum principle, i.e., it does not introduce new extrema. Note that the double contraction is applied as k v h γ,k . As commented above, we need to apply constraints to all hanging DOFs to keep conformity. The value of the hanging DOF needs to be equal to the value of the interpolation of the unknown at the neighboring coarser elements. That is, given i ∈ N hg h and its neighboring (coarse) In general, we will represent this constraint with v i = j∈N h C ij v j . For details of the implementation of this kind of constraints see [8][9][10]. Finally, to obtain the fully discrete problem, we consider a partition of the time domain (0, T ] into n ts sub-intervals of length (t n , t n+1 ]. Then, at every time step n = 0, . . . , n ts − 1, the discrete problem consists in solving The m × m-matrices relating nodes i, j ∈ N h are given by where Einstein summation applies, β, γ, ξ, η ∈ {1, . . . , m} are the component indices, δ βγ is the Kronecker delta, and M, K, and G are the terms arising from applying the constraints in the mass, flux and body forces terms.
2.3. Stability properties. Finally, let us review some concepts required for discussing the stabilization method used in the subsequent sections. Let us recall some definitions used for scalar problems.

Definition 2.3 (LED).
A scheme is local extremum diminishing if, for every u i that is a local discrete maximum (resp. minimum), is satisfied.
One possible strategy to satisfy the above properties consists in designing a scheme that yields a positive diagonal mass matrix and a stiffness matrix that satisfies In this case, it is possible to rewrite the system as As shown in [20,39], such a scheme satisfies the local DMP for steady problems and it is also local extremum diminishing (LED) when applied to transient problems. Stability properties for hyperbolic systems can be based on the extension of the above to hyperbolic systems in characteristic variables. In this direction, we define local bounds preserving schemes as follows.
is said to be local bounds preserving if M is diagonal with positive entries (i.e., M ij = m i δ ij I m×m ), A ij has non-positive eigenvalues for every j = i, and j A ij = 0.
Unfortunately, to the best of our knowledge, satisfying this definition does not ensure positivity of density, internal energy, or non-decreasing entropy. In any case, numerical schemes based on this definition have shown good numerical behavior [34,37,42,47].
Several stabilization strategies have been defined based on the previous ideas. One of the most simple strategies consists in adding a scalar artificial diffusion term proportional to the spectral radius of A ij [34,45]. This strategy is usually called Rusanov artificial diffusion, since the scheme results in the Rusanov Riemann solver for linear FEs in one dimension [34,56]. Without any special treatment, the resulting scheme is only first order accurate. The key for recovering high-order convergence is to modulate the action of the artificial diffusion term, and restrict its action to the vicinity of discontinuities. In the present work, our stabilization term for systems of equations is based on Rusanov artificial diffusion and a differentiable shock detector recently developed for scalar problems [3,15].
Finally, it is also important to define the concept of linearity preservation.
Definition 2.5. Given u h ∈ V h , a stabilization scheme is said to be linearity preserving if at any region such that u β h ∈ P 1 (Ω) ∀ β ∈ C, then B(u h ) = 0.

Nonlinear stabilization
In this section, we describe the additional terms used for the stabilized problem (2). In particular, we use the stabilization terms defined in [3] for the scalar problem and [5] for Euler: Here, is the graph-Laplacian operator defined as (i, j) . = 2δ ij − 1 (see [3,26]). In the case of a scalar problem, m = 1, the nodal artificial diffusion ν ij (w h ) is defined as We denote by α(w h ) the scalar shock detector used for computing the artificial diffusion parameter.
In the case of the Euler equations, the element-wise artificial diffusion ν e ij (w h ) is defined as where λ max ij is the spectral radius of the elemental convection matrix relating nodes i, j ∈ N h , i.e., ρ f (u ij ) · (∇ϕ j , ϕ i ) Ke , where u ij is the Roe average between u i and u j (see (6)). Notice that due to the usage of the FE approximation of the flux, λ max ik is the spectral radius of ρ f (u ij ) · (∇ϕ k , ϕ i ) Ke , where k ∈ N hg h . As previously introduced, this artificial diffusion term is based on Rusanov scalar diffusion [36]. It is important to mention that the eigenvalues of these matrices can be easily computed as where the velocity v ij and sound speed c are computed using the Roe averaged values We denote by α i (w h ) the shock detector used for modulating the action of the artificial diffusion term. The idea behind the definition of this detector is to minimize the amount of artificial diffusion introduced while stabilizing any oscillatory behavior. We ensure that Def. 2.4 is satisfied in regions where the local DMP is violated (see Def. 2.2) for any set of components. α i (w h ) must be a positive real number that takes value 1 when u h (x i ) is an inadmissible value of u h , and smaller than 1 otherwise. To this end, we define where C is the set of components that are used to detect inadmissible values of u h , e.g. density and total energy in the case of Euler equations. For simplicity, we restrict ourselves to the components of u h . However, derived quantities such as the pressure or internal energy can be also used. For scalar equations, since the stabilization is defined for the assembled system, the shock detector α i only needs to be defined for i ∈ N h . However, for the elemental definition used for Euler equations, it is also required for i ∈ N hg h . In that case, we use the maximum of its constraining nodes, i.e., Let us recall some useful notation from [3] to introduce the scalar shock detector α i (w h ). Let r ij = x j − x i be the vector pointing from node x i to x j with i, j ∈ N h andr ij . = rij |rij | . Recall that the set of points x j for j ∈ N h (Ω i )\{i} define the macroelement Ω i around node x i . Let x sym ij be the point at the intersection between ∂Ω i and the line that passes through x i and x j that is not x j (see Fig. 2 in two dimensions, let as call a and b the indices of the vertices such that they define the edge in ∂Ω i that contains x sym ij . We define u sym j as the value of u h at x sym ij , i.e., u h (x sym ij ).
In the present work, for each component in C, we use the same shock detector developed in [3]. Let us recall its definition We know from [3, Lm. 3.1] that (8) gets values between 0 and 1, and it is only equal to one if u β h (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 ∈ P 1 , the shock detector vanishes when the solution is linear. Thus, it is also linearly preserving for every component in C. This result follows directly from [3,Th. 4.5].
The final stabilized problem in matrix form reads as follows.
for n = 1, ..., n ts , where is the stabilization matrix. Let us show that adapted non-conforming meshes do not jeopardize any of the stability properties defined in Sect. 2.3 and proved for conforming meshes in [4,8].
Corollary 3.1 (DMP). The solution of the discrete problem (9) with m = 1 and using the shock detector (8) satisfies the local DMP in Def. 2.2 if g = 0 and, for every control point i ∈ N h such that u i is a local discrete extremum, it holds: Moreover, the resulting scheme is linearity-preserving as defined in Def. 2.5, i.e., B ij (u h ) = 0 for u h ∈ P 1 (Ω i ).
Proof. The stabilization scheme for scalar problems is defined on the assembled system. Hence, the modifications introduced in the assembly procedure do not affect the reasoning in the proof of [4, Thm.  Proof. If component β ∈ C of u h has an extremum at x i , we know from [3, Lm. 3.1] that α i (u β h ) = 1. Moreover, it is easy to see from (7) that α i (u h ) = 1. In this case, M ij (u h ) = δ ij j M ij . Hence, M ij (u h ) = 0 for j = i and M ii (u h ) = m i . Therefore, we can rewrite the system as follows We need to prove that the eigenvalues of K ij (u ij ) are non-positive. To this end, let us show that the following inequality holds From (5), it is easy to check that ρ f (u ij ) · (∇ϕ j , ϕ i ) Ke = |v ij · c e ij | + c ij c e ij . We have that where c e ij = (∇ϕ j , ϕ i ) Ke . Thus, Therefore, e ρ(K e ij (u ij )) ≥ ρ(K ij (u ij )). Moreover, by definition (see (4)), Furthermore, it is easy to infer from (3) that ρ(B e ij (u ij )) ≥ ρ(K e ij (u ij )). Hence, ρ(B ij (u ij )) ≥ ρ(K ij (u ij )). Finally, since K ij = K ij + B ij and B ij = e B e ij = e −ν e ij I m×m for all j = i. Then, the maximum eigenvalue of K ij (u ij ) is non-positive, which completes the proof.
Notice that it is essential to apply properly the constraints at the flux FE approximation, i.e., . Otherwise, it is not possible to formally prove local bound preservation. However, experimental results in the present work show that using f (u k ) does not affect the overall performance of the scheme.
3.1. Differentiable stabilization. In the case of steady, or implicit time integration, differentiability plays a role in the convergence behavior of the nonlinear solver. This is especially important if one wants to use Newton's method. The authors show in [3][4][5] that nonlinear convergence can be improved after few modifications to make the scheme twice-differentiable. In this section, we introduce a set of regularizations applied to all non-differentiable functions present in the stabilized scheme introduced above. In order to regularize these functions, we follow the same strategy as in [3][4][5]. Absolute values are replaced by Note that |x| 2,ε h ≤ |x| ≤ |x| 1,ε h . Next, we also use the smooth maximum function In addition, we need a smooth function to limit the value of any given quantity to one. To this end, we use The set of twice-differentiable functions defined above allows us to redefine the stabilization term introduced in Sect. 3. In particular, we definẽ whereν Let us note that λ max ij needs to be regularized as λ max The shock detector is also regularized as follows: In the case of the component shock detector we recall the definition in [3, Eq. 18] where ζ h is a small value for preventing division by zero. Finally, the twice-differentiable stabilized scheme reads: find u h ∈ V h such that u h = u h on ∂Ω, u h = u 0h at t = 0, and whereM Proof. For an extreme value of u β h , since |x| 2,ε h ≤ |x| ≤ |x| 1,ε h the quotient of (11) is larger than one. Hence, by definition of Z(x), α ε h ,i is equal to 1. At this point, it is easy to check thatν e ij ≥ ν e ij in virtue of the definition of max σ h . Therefore, ρ(B e ij (u h )) ≥ ρ(B e ij (u h )), completing the proof.
Moreover, it is important to mention that the differentiable shock detector is weakly linearlypreserving as ζ h tends to zero. This result follows directly from [3]. In order to obtain a differentiable operator, we have added a set of regularizations that rely on different parameters, e.g., σ h , ε h , ζ h . Giving a proper scaling of these parameters is essential to recover theoretic convergence rates. In particular, we use the following relations where for the scalar problem λ max is simply v , d is the spatial dimension of the problem, L is a characteristic length.

Adaptive mesh refinement
The motivation of an adaptive FE method is to solve (9) up to a certain tolerance (or resolution) using the minimum number of DOFs. To this end, the solution error (e h = u − u h ) is estimated at each element. With this information at hand, it is possible to iteratively adapt the resolution of the mesh at certain regions. This process can be divided into two parts: estimating the error at every cell, and deciding which and how many cells need to be refined or coarsened. This procedure is performed iteratively until a desired tolerance is achieved or, alternatively, a number of elements is reached. In the present work, we start with a rather coarse mesh and perform the following steps till reaching a stopping criterion: (1) Compute solution u h ; (2) Estimate the error e h ; (3) Select all cells that need to be refinement or coarsened; (4) Update the mesh, and project the solution to the new mesh.
In some cases, the refinement might be driven by features of the solution instead of a classical error estimator. For instance, one may decide to refine the regions around discontinuities. In this scenario, one could use a expression that does not estimate the error, but it allows to concentrate the elements around discontinuities.

Error estimators.
One of the keys of AMR is the ability to provide a good estimation of the error. Several error estimators have been proposed to date [1,22,29,30,58]. These can be classified, at least, in two main types. Some authors [22,29,49,50,52] try to compute an upper bound of the error for every cell. Then, provided a user defined tolerance, one can decide to refine or coarsen each cell. However, an adjoint problem needs to be solved in order to compute this upper bound [17,29]. It is possible to approximate the error bounds without solving an adjoint problem only for simple cases, see [29]. Therefore, this kind of error estimators increases the computational cost substantially. Alternatively, one can simply determine the distribution of the error in the mesh and use this information to drive an adaptivity algorithm. In this scenario, some authors [1,14,38,48,58,59] drive the adaptivity process with the solution gradient. In this case, explicit expressions of the estimated error are possible, requiring less computational resources than the previous option.
In general, the adaptive procedure can be described as follows. Given a finite element solution u h , the error e h is approximated as e h ≈ ∇u−∇u h . Then, the reconstruction is used as an approximation of the exact gradient. This strategy is based on superconvergence of special recovery techniques (see [59] and refs. in [1,48]). Kuzmin and co-workers [14,48] follow [58] to reconstruct an approximation of the exact gradient. Kelly et al. [30] proposed a well-known estimator based on gradient recovery: where η K is the estimated error at every element, K. The main advantages of this estimator are its simplicity and its low computational cost. For these reasons, this estimator is used in the present work. It is worth mentioning that our problems of interest are characterized by exhibiting discontinuities, where the error concentrates. These regions are susceptible to develop instabilities, and thus, these are the regions in which the shock capturing is activated. Therefore, it is natural to use the shock capturing to drive the adaptivity procedure. We propose an estimator based on the graph Laplacian (i, j) present in the stabilization term (3). This way, we reuse available information and reduce the computational overhead associated with error estimation. The indicator reads: where β ∈ C is the index of the specific component analyzed. This estimator is expected to yield high values around shocks and low values in smooth regions.

Refinement strategy.
After the error has been estimated for every element, one needs to decide which element needs to be refined and which one coarsened. If an upper bound of the error is computed, then one may use a given tolerance to make this decision. However, in the present case this is not available. A classical alternative is to refine/coarsen a fixed amount of elements at every iteration [9,11]. In the present study, a 30% of the elements with higher error estimates are refined whereas a 10% of the elements with lower error estimates are coarsened. This percentages are arbitrary and other choices are valid. Notice that using this setting in two dimensions the number of elements is almost doubled at every iteration. We make use of the parallel nth element algorithm [9,55] to efficiently determine the error estimator thresholds for refining or coarsening the elements.

Nonlinear solver
In this section, we describe the method used for solving the nonlinear system of equations arising from the scheme introduced above. In particular, we use a hybrid Picard-Newton approach in order to increase the robustness of the nonlinear solver. Moreover, we also make use of a line-search method to improve the nonlinear convergence.
We define the residual of the equation (12) at the k-th iteration as Hence, the Jacobian is defined as Therefore, Newton method consists in solving J(u k,n+1 h )∆U k+1,n+1 = −R(u k,n+1 h ). It is well known that Newton method can diverge if the initial guess of the solution u 0,n+1 h is not close enough to the solution. In order to improve robustness, we use a line-search method to update the solution at every time step. The new approximation is computed as U k+1,n+1 = U k,n+1 + λ∆U k+1,n+1 , where λ is obtained using a standard cubic backtracking algorithm.
As introduced at the beginning of the section, we also use a hybrid approach combining Newton method with Picard linearization. Picard nonlinear iterator can be obtained removing the last two terms of (13), i.e., Clearly, it is equivalent to Moreover, we modify the left hand side terms in (14); we use α i = 1 for computing these terms while we use the value obtained from (7) for the residual. Using this strategy, the solution remains unaltered but the obtained approximations u k,n+1 h for intermediate values of k are more diffusive. Even though this modification slows the nonlinear convergence, it is essential at the first iterations. Otherwise, the robustness of the method might be jeopardized.
The resulting iterative nonlinear solver consists in the following steps. We iterate Picard method in (14), with the modification described above, until the L 2 norm of the residual is smaller than a given tolerance. In the present work, we use a tolerance of 10 −2 . Afterwards, Newton method with the exact Jacobian in (13) is used until the desired nonlinear convergence criteria is satisfied. We summarize the nonlinear solver introduced above in Alg. 1.

Numerical results
In this section, we perform several numerical experiments to assess the numerical scheme introduced in the previous sections. First, we perform a convergence analysis to assess its implementation. Then, we use steady benchmark tests to analyze the effectiveness of the high-order scheme in the context of AMR. In particular, we compare the nonlinear scheme in (12) with its linear (first order) counterpart, i.e., using α ε h ,i (u h ) ≡ 1.
6.1. Convergence. First, the convergence to a discontinuous solution is analyzed. To this end, we solve two different problems. On the one hand, the following scalar problem is solved where v(x, y) . = ( 1 /2, sin −π /3), and inflow boundary conditions u D = 1 on {x = 0} ∩ {y > 0.7} and y = 1, while u D = 0 at the rest of the inflow boundary. This problem has the following analytical solution For the Euler equations, the problem is the well known compression corner test [2,36], also known as oblique shock test [51,53]. This benchmark consists in a supersonic flow impinging to a wall at an angle. We use a [0, 1] 2 domain with a M = 2 flow at 10 • with respect to the wall. This leads to two flow regions separated by an oblique shock at 29.3 • , see Fig. 3.
Since the solution is not smooth, we expect linear convergence rates in the L 1 -norm. Fig. 4 shows the convergence behavior of both problems with uniform mesh refinements. The experimental convergence rate measured for the scalar transport problem is 0.82, whereas the convergence rate measured is 0.94 for the compression corner test. Therefore, both tests exhibit the expected convergence behavior.
6.2. Linear discontinuity. For this test, we use again the problem in (15). The purpose of this test is twofold. On the one hand, we analyze the effectiveness of the proposed estimator. On the other hand, we compare the effectiveness of the linear and nonlinear stabilization methods. Specifically, this    effectiveness is measured as follows. For a given error, we consider a method more effective if it requires less computational time, independently of the number of elements required. In addition, we also solve the problem for successive uniformly refined meshes in order to evaluate the effect of AMR. For all comparisons, we start with a coarse mesh of 16 × 16 elements, and proceed adapting the mesh up to a maximum number of elements. For the nonlinear stabilization, we set a maximum of 10 4 elements. The maximum number of elements for the low-order method is 10 5 . The uniform mesh is refined up to a 256 × 256 mesh. We use a nonlinear tolerance of ∆u h / u h < 10 −4 , and a maximum of 500 iterations. Fig. 5 shows the evolution of the AMR algorithm for both error estimators. The results shown in this picture have been obtained using the linear stabilization, and the left-most column using the nonlinear one. It can be observed that both Kelly (η K ) and graph Laplacian (η K ) estimators refine in the vicinity of the shock. However, the graph Laplacian operator clearly outperforms Kelly estimator.
Figs. 6-8 compare the effectiveness of the low-order and the high-order stabilization schemes. The results are obtained for the stabilization parameter q = 1, q = 2, and q = 10, respectively.
At Fig. 7, the nonlinear stabilization is able to converge the nonlinear problem efficiently and the overhead of solving a nonlinear problem does not strongly affect the overall performance. We note that for the linear scheme the problem is linear. It can be observed that the convergence rate (against time) is much higher for the nonlinear scheme. The linear scheme requires less computational time for coarser meshes but the nonlinear scheme is more effective for tighter accuracies.
We can observe in Fig. 8 the convergence problems of the nonlinear stabilization at some steps of the refinement procedure. Even though using q = 10 improves the accuracy of the method, it also increases the computational cost since the nonlinear problem is harder to solve. This results in an efficiency slightly lower than the linear stabilization. Figure 5. Evolution of the mesh refinement process.η K with high-order scheme is used in the left column. Low-order scheme with Kelly estimator is used in the central column.η K with low-order scheme is used in the right column. For the low-order scheme from top to bottom results have been obtained at refinement step 1, 2, 3, 9, and 9. For the high-order refinement steps are 1, 2, 3, 5, and 5.       The analytical solution of this particular configuration consists in the transport of the inflow profile in the direction of the convection. As a result, the solution at the outflow boundary, corresponding to y = 0, is u(x, 0) = u(0, x). We start with a coarse mesh of 16 × 16 elements in all cases, and proceed adapting the mesh up to a maximum number of elements. For the nonlinear stabilization, we set a maximum of 5 · 10 3 elements. The maximum number of elements for the linear stabilization is 5 · 10 4 . We use a nonlinear tolerance of ∆u h / u h < 10 −4 , and a maximum of 500 iterations.        Low order, Kelly, η High order, Kelly, η Low order,η High order,η Figure 11. Time and elements convergence comparison for the transport problem with a circular convection field, q = 10. Figure 12. Evolution of the mesh refinement process.η K with high-order scheme is used in the left column. Low-order scheme with Kelly estimator is used in the central column.η K with low-order scheme is used in the right column. For the low-order scheme from top to bottom results have been obtained at refinement step 1, 2, 3, 7, and 7. For the high-order refinement steps are 1, 2, 3, 4, and 4. Figs. 9-11 compare the effectiveness of the linear and nonlinear stabilization. These results use the stabilization parameter q = 1, q = 2, and q = 10, respectively. In Fig. 10, the high-order scheme is able to converge efficiently and the overhead of solving a nonlinear problem does not strongly affect the overall performance. Nevertheless, the low-order scheme requires less computational time for any given error. However, it can be observed that the convergence rate (in time) is much higher for the high-order scheme. Therefore, it is expected to outperform the low-order scheme for more refined meshes.
In contrast, we do not observe the significant convergence problems in Fig. 11 even though the linear stabilization is slightly more efficient. Again, the convergence rate (versus time) is higher for the nonlinear stabilization and it will be better for more refined meshes. Fig. 12 shows the evolution of the AMR algorithm for both η K andη K estimators with the linear stabilization andη K with the nonlinear one. It can be observed that both Kelly (η K ) and graph Laplacian indicators detect the regions that require more resolution. In any case, as in the previous example, the graph Laplacian operator (η K ) performs slightly better.
6.4. Compression corner. Let us consider now the Euler equations. We start with the compression corner test (see Fig. 3). We analyze the effectiveness of the high-order scheme, and evaluate the performance of the graph Laplacian estimator. We start with a coarse mesh of 16 × 16 elements, and adapt it up to a maximum number of elements. For the high-order method, we set a maximum of 5 · 10 3 elements. The maximum number of elements for the low-order method is 5 · 10 4 . We use a nonlinear tolerance of ∆u h / u h < 10 −4 and a maximum of 500 iterations. Fig. 13 compares the effectiveness of the low-order and the high-order stabilization schemes for different values of q. The high-order scheme is able to converge efficiently and the overhead of solving a nonlinear problem does not affect the overall performance. In this case, the low-order and the highorder schemes require similar computational time for any given error. Actually, for the finer meshes, the high-order scheme with either q = 1 or q = 2 already performs better than the low-order scheme. However, for some meshes the high-order scheme exhibits convergence problems. In the case of q = 10 the cost of converging the nonlinear problem does not compensate the increase in computational cost.
In Fig. 14, we depict the refinement evolution for the graph Laplacian estimator (η K ) for linear and nonlinear stabilization. As expected, we can observe that for the high-order method the scheme is able to resolve the shock with less refinement steps. The linear stabilization is able to provide well-resolved shocks at the final refinement step.  Low order,η High order,η, q = 1 High order,η, q = 2 High order,η, q = 10  Figure 14. Evolution of the mesh refinement process.η K with high-order (right) and low-order (left) schemes are used. For the low-order scheme from top to bottom results have been obtained at refinement step 1, 2, 3, 8, and 8. For the high-order refinement steps are 1, 2, 3, 4, and 4. leads to a steady shock separating both flow regimes that is reflected at the wall producing a third different flow state behind it. A sketch of this benchmark test is given in Fig. 15. The flow states at each region have been collected in Tab. 1.  We analyze the effectiveness of the high-order scheme, and evaluate the performance of the graph Laplacian estimator. We start with a coarse mesh of 16 × 64 elements and adapt the mesh till a certain number of elements is reached. For the high-order method, we set a maximum of 10 4 elements. The maximum number of elements for the low-order method is 3 · 10 5 . We use a nonlinear tolerance of ∆u h / u h < 10 −4 and a maximum of 500 iterations. Fig. 16 compares the effectiveness of the low-order and the high-order stabilization schemes for different values of q. The high-order scheme converges efficiently and the overhead of solving a nonlinear problem does not affect the overall performance. Actually, for the most refined meshes the high-order method is more efficient than the low-order one. As for the previous problem, Fig. 16 shows that the high-order scheme can present nonlinear convergence problems at some steps of the refinement process. However, as the mesh becomes more adapted to the problem this issue is reduced.  Low order,η High order,η, q = 1 High order,η, q = 2 High order,η, q = 10 Figure 16. Time and elements convergence comparison for the reflected shock problem.
In Fig. 17 we depict the refinement evolution for the graph Laplacian estimator (η K ) for the low-order scheme. In these figures it can be observed how the graph Laplacian estimator is able to concentrate all the resolution at the shock location. Finally, we can conclude from the lower two figures that both schemes resolve the shocks properly after the mesh has been refined enough. Figure 17. Evolution of the mesh refinement process.η K with low-order scheme is used. For the low-order scheme from top to bottom results have been obtained at refinement step 1, 2, 3, 4, 5, 6, and 7. The lower two figures are the high-order (top) and low-order (bottom) results at their last refinement step.

Conclusions
The stabilization schemes in [3,5] have been extended and assessed in the AMR context for nonconforming hierarchical octree meshes. The work focuses in assessing the effectiveness of linear (first-order) and nonlinear (higher-order) stabilization. We focus the comparison in terms of accuracy versus computational time.
The results indicate that linear stabilization is more effective for coarse meshes. In this case, the computational cost required to solve the stiff nonlinear problem due to the nonlinear stabilization does not compensate the improvement in the accuracy. This is especially evident for linear systems of partial differential equations (PDEs). On the contrary, as the mesh is refined and properly adapted to the shocks, nonlinear stabilization pays the price. Even though increasing the value of q in the nonlinear stabilization (a parameter that makes shocks sharper but hinders nonlinear convergence) improves accuracy, it turns to be far more effective to refine the mesh further for low values of q. Nevertheless, it is worth mentioning that high-order method might exhibit nonlinear convergence problems for some meshes.
In addition, a new error estimator have been proposed. The proposed estimator is based on the graph Laplacian used in the definition of the stabilization method. Numerical results show that this shock detector is able to perform better that the well known Kelly estimator for problems with shocks or discontinuities.