Compressible flows can be characterised by abrupt changes in density within the flow domain often referred to as shocks. These discontinuities can lead to numerical instabilities (Gibbs phenomena). This problem is prevented by locally adding a diffusion term to the equations to damp the numerical oscillations.

For the non-smooth artificial viscosity model the added artificial viscosity is constant in each element and discontinuous between the elements. The Euler system is augmented by an added Laplacian term on right hand side of equation 9.1 [35]. The diffusivity of the system is controlled by a variable viscosity coefficient ε. For consistency ε is proportional to the element size and inversely proportional to the polynomial order. Finally, from physical considerations ε needs to be proportional to the maximum characteristic speed of the problem. The final form of the artificial viscosity is

| (9.8) |

where S is a sensor.

As shock sensor, a modal resolution-based indicator is used

| (9.9) |

where ⟨⋅,⋅⟩ represents a L^{2} inner product, q and are the full and truncated expansions of a
state variable (in our case density)

| (9.10) |

then the constant element-wise sensor is computed as follows

| (9.11) |

where s_{0} = s_{κ} - 4.25 log_{10}(p).

To enable the non-smooth viscosity model, the following line has to be added to the
`SOLVERINFO`

section:

The diffusivity and the sensor can be controlled by the following parameters:

A sensor based p-adaptive algorithm is implemented to optimise the computational cost and
accuracy. The DG scheme allows one to use different polynomial orders since the fluxes over
the elements are determined using a Riemann solver and there is now further coupling
between the elements. Furthermore, the initial p-adaptive algorithm uses the same
sensor as the shock capturing algorithm to identify the smoothness of the local
solution so it rather straightforward to implement both algorithms at the same
time.

The polynomial order in each element can be adjusted based on the sensor value that is
obtained. Initially, a converged solution is obtained after which the sensor in each element is
calculated. Based on the determined sensor value and the pre-defined sensor thresholds, it is
decided to increase, decrease or maintain the degree of the polynomial approximation in each
element and a new converged solution is obtained.

| (9.12) |

For now, the threshold values s_{e}, s_{ds}, s_{sm} and s_{fl} are determined empirically by looking at the
sensor distribution in the domain. Once these values are set, two .txt files are outputted, one
that has the composites called VariablePComposites.txt and one with the expansions called
VariablePExpansions.txt. These values have to copied into a new .xml file to create the
adapted mesh.

Aliasing effects, arising as a consequence of the nonlinearity of the underlying problem, need to be address to stabilise the simulations. Aliasing appears when nonlinear quantities are calculated at an insufficient number of quadrature points. We can identify two types of nonlinearities:

- PDE nonlinearities, related to the nonlinear and quasi-linear fluxes.
- Geometrical nonlinearities, related to the deformed/curves meshes.

We consider two de-aliasing strategies based on the concept of consistent integration:

- Local dealiasing: It only targets the PDE-aliasing sources, applying a consistent integration of them locally.
- Global dealiasing: It targets both the PDE and the geometrical-aliasing sources. It requires a richer quadrature order to consistently integrate the nonlinear fluxes, the geometric factors, the mass matrix and the boundary term.

Since Nektar++ tackles separately the PDE and geometric aliasing during the projection and solution of the equations, to consistently integrate all the nonlinearities in the compressible NavierStokes equations, the quadrature points should be selected based on the maximum order of the nonlinearities:

| (9.13) |

where Q_{min} is the minimum required number of quadrature points to exactly integrate the
highest-degree of nonlinearity, P_{exp} being the order of the polynomial expansion and P_{geom}
being the geometric order of the mesh. Bear in mind that we are using a discontinuous
discretisation, meaning that aliasing effect are not fully controlled, since the boundary terms
introduce non-polynomial functions into the problem.

To enable the global de-aliasing technique, modify the number of quadrature points by:

1<E COMPOSITE="[101]"

2 BASISTYPE="Modified_A,Modified_A"

3 NUMMODES="7,7"

4 POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre"

5 NUMPOINTS="14,14"

6 FIELDS="rho,rhou,rhov,E"

7/>

2 BASISTYPE="Modified_A,Modified_A"

3 NUMMODES="7,7"

4 POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre"

5 NUMPOINTS="14,14"

6 FIELDS="rho,rhou,rhov,E"

7/>

where `NUMMODES`

corresponds to P+1, where P is the order of the polynomial used to
approximate the solution. `NUMPOINTS`

specifies the number of quadrature points.