$rp$-adaptation for compressible flows

We present an $rp$-adaptation strategy for the high-fidelity simulation of compressible inviscid flows with shocks. The mesh resolution in regions of flow discontinuities is increased by using a variational optimiser to $r$-adapt the mesh and cluster degrees of freedom there. In regions of smooth flow, we locally increase or decrease the local resolution through increasing or decreasing the polynomial order of the elements. This dual approach allows us to take advantage of the strengths of both methods for best computational performance, thereby reducing the overall cost of the simulation. The adaptation workflow uses a sensor for both discontinuities and smooth regions that is cheap to calculate, but the framework is general and could be used in conjunction with other feature-based sensors or error estimators. We demonstrate this proof-of-concept using two geometries at transonic and supersonic flow regimes. The method was implemented in the open-source spectral/$hp$ element framework Nektar++, and its dedicated high-order mesh generation tool NekMesh. The results show that the proposed $rp$-adaptation methodology is a reasonably cost-effective way of improving accuracy.


INTRODUCTION
The accurate and high-fidelity simulation of high-speed compressible flows is, at present, a problem of significant interest to the aeronautics community, particularly in relation to aviation in which such conditions are routinely encountered. The complex and interdependent fluid phenomena found in this regime pose a difficult challenge for numerical modelling, with a stark contrast between regions of smooth flow, boundary layers near solid walls where large velocity gradients are present, and the interaction with shock waves or shear layers where the fluid properties change sharply in a discontinuous manner.
The use of high-order spectral/ℎ element methods in the simulation of compressible fluid dynamics is now becoming increasingly common for high-fidelity large-eddy simulations and direct numerical simulations of realistic aeronautical configurations (de Grazia, Moxey, Sherwin, Kravtsova, & Ruban 2018). As in traditional low-order methods, the domain of interest is partitioned into finite elements; however, these elements are also equipped high-order polynomial expansions, as opposed to traditional linear shape functions. This yields several advantages in terms of computational performance, as well as enhanced numerical resolution as is increased. However, in the presence of shocks and discontinuities, the latter advantage will not be realised, and can lead to significant issues in terms of stability and accuracy in the resolution of shocks.
A common approach used in the resolution of discontinuous features is to refine these regions in an adaptive manner, so that the mesh resolution around the features is increased. In broad terms, the error of a computed solution which is sufficiently 0 Abbreviations: BC: boundary conditions; DOF: degrees of freedom; HPC: high-performance computing; I/O: input-output arXiv:1909.10973v1 [physics.comp-ph] 24 Sep 2019 smooth can be roughly expressed as ≈ ℎ , where is a constant related to the measure of the solution regularity, ℎ is the mesh size, and is the polynomial order. Mesh adaptation is concerned with achieving increased resolution by either locally reducing the mesh size, ℎ, or locally increasing the polynomial order, . Due to its higher convergence rates, p-adaptation is typically preferred over h-adaptation for smooth flow regions as discussed by Burbeau and Sagaut (2005), Li, Premasuthan, and Jameson (2010) and Ekelschot, Moxey, Sherwin, and Peiró (2016), whereas the opposite is true where flow discontinuities exist. The reason for the latterh-adaptation being preferred for flow discontinuities -is that the representation of shocks by highorder discretisations leads to numerical oscillations that must be smoothed out by the addition of high-order dissipation terms. This effectively means that the high-order degrees of freedom (DOF) are wasted in the vicinity of shocks.
To address these issues, in this article we present a proof-of-concept strategy based on rp-adaptation to best take advantage of h-type, through r-adaptation, and p-type local resolution modifications. The r-adaptation procedure is based on the work of Marcon, Turner, Moxey, Sherwin, and Peiró (2017) where a variational optimiser is used to deform the mesh. By targeting a small element size in regions of shocks, the optimiser deforms the mesh and clusters nodes in said regions. By effectively redistributing DOF, h-type refinement is obtained at flow discontinuities. We then apply p-adaptation to this adapted mesh, based on the work of Ekelschot et al. (2016), to better resolve regions of smooth flow. Throughout the work, we focus on the simulaton of inviscid flows and focus on the challenge of efficiently modelling smooth flows with embedded discontinuities.
In both procedures, we use a discontinuity sensor introduced by Persson and Peraire (2006). This sensor is easily computed and essentially looks at the energy of the higher modes to determine the level of resolution of the solution. The purposes of this sensor are three-fold: first, it adds artificial viscosity to the governing equations, based on values of the sensor, to stabilise the solution in the presence of shocks; second, it identifies regions of flow discontinuities based on values of the sensor, as used for the artificial viscosity, to drive r-adaptation; and third, it locally increases or decreases the local polynomial approximation based on the values of the sensor.
We present this proof-of-concept methodology as follows. Sect. 2 introduces the governing equations in continuous and discrete forms. Sect. 3 describes the spectral/hp discontinuous Galerkin discretisation used in Nektar++, with Sect. 3.1 covering the formulation of the discontinuity sensor and the artificial viscosity. Sect. 4 recalls the work of Marcon et al. (2017) on variational r-adaptation. Sect. 5 summarises the p-adaptation strategy originally proposed by Ekelschot et al. (2016). Sect. 6 describes the novel dual rp-adaptation workflow. Finally, we present two numerical examples in Sect. 7: a transonic flow past a NACA 0012 profile at a free-stream Mach number of 0.8 an incidence of 1.25 degrees, and a supersonic flow at a free-stream Mach number of 3 past an engine intake that exhibits a complex shock pattern in its diffuser.

GOVERNING EQUATIONS
The Euler equations of inviscid compressible flow are written, in a two-dimensional Cartesian frame of reference with coordinates = ( 1 , 2 ) within a domain Ω with boundary Γ, as Here = [ , 1 , 2 , ] is the vector of conserved variables, where is the density, the Cartesian components of the velocity are = ( 1 , 2 ), and is the total energy. The terms and denote the convective and dissipative fluxes, respectively. A dissipative flux is required to stabilise the solution in the presence of shocks which is chosen to be of the form where is an artificial viscosity coefficient that will be discussed in detail in Sect. 3.1. The components of the convective flux, = ( 1 , 2 ), are given by where is the total enthalpy and is the pressure. The total enthalpy is defined as and, to close the system, the pressure for a perfect gas is given by where is the ratio of specific heats and its value for air is = 1.4. The setting of the problem is completed through a suitable choice of initial and boundary conditions. Given that only steadystate problems are of interest here, we start the simulation with a uniform flow at the given free-stream Mach number and flow incidence. Solid walls are modeled through the no-flow condition, ⋅ = 0, where denotes the wall outer normal. Far-field boundaries are weakly imposed through the normal boundary fluxes by specifying free-stream conditions, = ∞ , outside the boundary and evaluating the normal fluxes through a Riemann solver that accounts for the propagation of information across the boundary.

DISCONTINUOUS GALERKIN DISCRETISATION
To obtain a discrete solution of (1) via a high-order spectral/ℎ discontinuous Galerkin discretisation, we assume that the computational domain, Ω, is subdivided into non-overlapping elements, so that Ω = ⋃ =1 Ω , where ⋂ =1 Ω = ∅. We adopt a mixed formulation similar to that proposed by Bassi and Rebay (1996) and write (1) as We seek a discrete approximation within an element, Ω , of the form where ( ); = 1, … , , represent the elemental expansion functions for the high-order spectral/ℎ discontinuous Galerkin method available in Nektar++ and presented by Cantwell et al. (2015) and Moxey et al. (2019). Both the solution and test functions are discontinuous at the interface between elements. Following the standard Galerkin procedure, a weak form of the mixed formulation (6)- (7) is obtained as follows. The discrete version of equation (6) reads Using an approximation of the form (8) for both ℎ and ℎ , and applying Gauss' theorem this equation becomes where Γ denotes the boundary faces of the element Ω . The solution of this equation give us the discrete values of the first-order derivatives ℎ . To evaluate the integral expressions, we use an auxiliary mapping ∶  → to transform the local element coordinates = ( 1 , 2 ) to reference element coordinates = ( 1 , 2 ) such that −1 ≤ 1 , 2 ≤ 1 and all required operations take place in the reference element Ω st , see Fig. 1.
The weak form of equation (7) is obtained in a similar fashion to give The solution is discontinuous at the interface between the elements and the integrand in the boundary integral of (11) is substituted by a numerical flux function. The convective normal flux at an interface is approximated by a numerical flux calculated via a Riemann solver ⋅ Γ ≈  ( , + ; ) where + and are the values of the conservative variables on the external and internal sides of the interface with respect to the -th element. This mechanism allows information to pass from one element to the other. The evaluation of the diffusive normal flux at the interface follows the local discontinuous Galerkin (LDG) formulation proposed by Cockburn and Shu (1998), where it is approximated by and similarly where 12 = 1 2 , and 11 is an order 1 constant. The average and jump operators are defined as

Shock capturing via a discontinuity sensor
Our DG discretization of the Euler equations requires the addition of the diffusion flux, to stabilize the solutions in the presence of shock waves. The term in (2) is an artificial viscosity coefficient that allows dissipation to be selectively applied to shocks. For consistency ∼ ℎ∕ , Barter and Darmofal (2010) suggest where max = | | + is the local maximum wave speed of the system. The characteristic cell length ℎ is chosen as the minimum edge length of an element. Finally, for the artificial viscosity to vanish outside shocks it needs to be proportional to a shock sensor, , such as where 0 = (1) is a constant. To build the shock sensor, we adopt the modal resolution-based indicator proposed by Persson and Peraire (2006) which is element-wise constant and defined via an intermediary term where ⟨⋅, ⋅⟩ represents a 2 inner product, and̃ are the full and truncated expansions of a state variable where are the basis functions,̂ the associated coefficients, and ( ) the size of the expansion of order . In our case the test variable is chosen to be the density . To spatially smooth out the variation of the values of the sensor, the constant element-wise sensor is computed as follows with 0 ∼ log 10 ( 4 ) from an analogy to Fourier coefficients decaying as 1∕ 2 , and needs to be sufficiently large to obtain a smooth shock profile. We select 0 = − − 4.25 log 10 ( ) , where and can be adjusted for a specific problem. We describe how to choose these parameters in Sect. 6.

4 R-ADAPTATION
In r-adaptation we are aiming at increasing resolution by locally reducing the mesh size, ℎ, whilst keeping the number of DOF in the mesh constant. This effectively requires us to cluster mesh nodes in the vicinity of those regions where additional resolution is required, e.g. shocks. We propose to accomplish this by adapting a variational framework for the optimisation of high-order meshes proposed by Turner, Peiró, and Moxey (2017).

Variational mesh optimisation
The objective of the variational framework by Turner et al. (2017) is to improve the quality of high-order curvilinear elements using a node-based optimisation approach using a formulation based on the energy of deformation. An important aspect of such energy-based formulation is that a suitable choice of the energy functional, namely polyconvex, would guarantee the existence of a minimum and therefore of a solution to the minimisation problem. Fig. 1 shows that a mapping exists from a reference element to a curvilinear high-order element . We can further decompose the mapping into two distinct mappings: a mapping from reference to ideal elements and a mapping from the ideal to the curvilinear elements. We define the ideal element as the high-order linear element, which after minimisation will be the element that the optimiser seeks to achieve.
Existing mappings between the reference, the ideal and the curvilinear elements. The ideal and curvilinear elements become respectively the target and adapted elements in the framework of r-adaptation.
From this ideal element, we calculate the deformation energy. The mesh is deformed to minimise an energy functional (∇ ): where (∇ ) is a formulation of the deformation energy. Turner et al. (2017) tested several formulations and found that the best results were obtained when using a hyperelastic model. For this model, the strain energy takes the form where and are material constants, is the right Cauchy-Green tensor, 1 is its trace and is the determinant of the Jacobian matrix = ∇ . We use this formulation in the work that follows. When optimising a high-order mesh, the ideal element of Fig. 1 corresponds to the linear element before generation of a high-order element by addition and projection of high-order nodes. To obtain best performance, not only high-order nodes but vertices are optimised too. Marcon et al. (2017) have proposed to change this ideal element and make it an arbitrary target element. The optimiser now aims at adapting element towards a size and a shape similar to the target element .

Improving the shock resolution via element scaling
Although in principle the target element can take any shape and size, we have adopted a practical approach in this work that aims at avoiding too large deformations. The rationale for this is that the definition of a target element that is very different from the ideal element -i.e. the initial linear element before r-adaptation -introduces extra energy in the system that the optimiser has to minimise and thus slows down the process. For this reason, we define a target element with respect to the ideal element . Marcon et al. (2017) demonstrated that the ideal element can be manipulated anisotropically by applying a metric tensor to the Jacobian of the mapping, = ∇ . We transform the ideal element into the target element through = We do not consider directionality in this work and only shrink elements where additional resolution is required, i.e. in the shock regions. In this case, the Jacobian is simply scaled by a linear shrinking factor . The metric tensor is simplified to to obtain This framework was implemented in NekMesh, an open-source software solution for the generation of geometry-accurate high-order meshes, part of the Nektar++ platform presented by Cantwell et al. (2015) and Moxey et al. (2019).

P-ADAPTATION
To enhance resolution in regions of smooth flow through local p-adaptation, we follow the procedure laid out by Ekelschot et al. (2016) and . The strategy is fairly straightforward. We increase the local resolution by increasing the polynomial order within the elements where the local error is estimated to be high and we decrease the local resolution or, equivalently, the elemental polynomial order within those elements where the local error is estimated to be low. We summarise this procedure in Algorithm 1 where denotes an individual element, and are its associated error indicator and polynomial order, and are the upper and lower error thresholds, and max and min are the maximum and minimum polynomial orders allowed. In this work, we use the formulation of sensor in equation (19) as the indicator of the discretisation error.

WORKFLOW FOR RP-ADAPTATION
Finally we attempt to combine the best properties of the two previous strategies in the same simulation to maximize their effect in increasing the resolution of compressible flow simulations. More specifically, r-adaptation will be responsible for the resolution of shocks whereas p-adaptation will resolve smooth flow regions. In the proposed rp-adaptation workflow these adaptative techniques will be alternatively applied in a sequence of steps that is described in the following.
We first generate an initial high-order mesh for the domain. We anticipate the requirements of r-adaptation and the need for DOF to be moved around when deforming the mesh. For this reason, we generate a relatively coarse mesh, but with enough resolution to allow for nodes movement. We then proceed to run the solver on this initial mesh and obtain a flow solution which represents our base solution. During this step, we have to determine appropriate parameter values for the artificial viscosity. As is common practice in codes based on artificial viscosity, we start with Nektar++ default values, and then tune the parameters for our specific problem. We adjust the artificial viscosity parameters ( and ) to ensure that artificial viscosity is only triggered in the direct vicinity of shock waves, and adjust the level of artificial viscosity ( 0 ) so that the shock is stable but not overly dissipative.
From this base solution, we apply r-adaptation to the mesh. We first extract the list of elements where artificial viscosity was added during the initial simulation. If the run was set up properly, these elements only represent the regions where a shock is present. From these elements, we extract their barycentres and assign an isotropic shrinking factor to them (see Sect. 4.2). For all the other elements, we also extract the barycentres and assign a factor = 1. In practical terms, we force elements in the shock regions to shrink and pull mesh nodes from all other parts of the mesh. This field of r factors is then supplied to the variational r-adaptation code which is then run on the linear mesh. The variational framework optimises the mesh so that each element is as close as possible to its target size, effectively moving nodes from areas of = 1 to areas of < 1. We also note that r-adaptation is run on the linear mesh before making the adapted mesh high-order again. This significantly speeds up the optimisation procedure and improves the validity of the final mesh. We then run the solver on the adapted high-order mesh and obtain a new solution with enhanced shock resolution. This procedure can optionally be repeated based on the new solution.
From this solution on the adapted mesh, we can run p-adaptation as described in Sect 5. At the end of each cycle, a sensor value is computed for each element and the local polynomial order of that element is decreased, kept the same or increased based on the value of the sensor. In this work, we use Nektar++ default values for the adaptation parameters: = −6, max = 6, = −8 and min = 2. The simulation then proceeds onto a new cycle and the process is repeated until a steady solution is obtained and the local polynomial orders do not vary.
We also allow the user a choice to restrict the polynomial order of elements within the shock regions. These are the zones that have been previously identified in the r-adaptation procedure. The local polynomial order of the elements in these regions is then set to a user-defined value, which should be typically low ( < 3). The proposed rp-adaptation workflow is summarised in Algorithm 2.

Algorithm 2
The rp-adaptive workflow.
generate an initial high-order mesh calculate the steady-state solution repeat extract shock areas based on sensor values apply r-adaptation in shock areas to linear mesh and re-project to high-order calculate the steady-state solution until shocks are well captured apply p-adaptation as described in Algorithm 1 calculate the final solution

NUMERICAL EXAMPLES
In this section, we use two different test cases to demonstrate the rp-adaptation workflow: a NACA 0012 aerofoil in transonic regime and a supersonic intake. Different difficulties arise for each of these test cases as we will discuss below. Most importantly, we take slightly different approaches when it comes to r-adaptation and p-adaptation.
For the first test, only two simple shocks are present on the upper and lower surfaces of the aerofoil. It is easy for the variational framework to pull mesh nodes from the smooth regions and we therefore run a single round of r-adaptation. For the second test, a complex diamond-like pattern of oblique shocks is expected due to reflections inside the internal configuration of the intake. The number of nodes available for deforming the mesh is limited, which places additional stress on the variational optimisation process. We therefore decide to take a two-step approach to r-adaptation where each step is run with a milder shrinking factor in order to retain good mesh quality.
In terms of p-adaptation, we use the first test case as a benchmark for our order restriction strategy. We study three different approaches to order restriction. In the first, no restriction is applied and the local polynomial order of elements in shock areas is left free to increase. In the second, we keep the local polynomial order of these elements constant, i.e. the order of the initial simulation. Finally, in the third case we immediately decrease the local polynomial order of these elements to the minimum allowed in the run.

NACA 0012
We first demonstrate the new technology on a canonical aeronautical test case: a transonic ( = 0.8) inviscid flow past a NACA 0012 aerofoil at 1.25 • angle of attack. Vassberg and Jameson (2010) have shown that this configuration produces two shocks: a strong shock on the suction side and a weak shock on the pressure side at approximately 60% and 35% of the chord respectively. This provides a relatively easy test case to showcase the technology where the shocks are quasi-vertical. The main difficulty lies in the relative weakness of the shock on the pressure side and in capturing it appropriately.
The domain used has external boundaries at a distance of 40 from the aerofoil, where denotes the chord length. We discretise the domain uniformly along the chord with an element size of ∼ 0.5 on the aerofoil boundary and create a smooth progression towards an element size of ∼ 10 on the outer boundary. The automatic sizing of elements in the field is determined through an octree system as described by Turner, Moxey, Sherwin, and Peiró (2016). The mesh is curvilinear of order = 4 and it is optimised using the variational framework described in Sect. 4.1. Fig. 2a (left) shows what the mesh looks like in the near field. The starting mesh is relatively coarse but it is run through the solver at uniform = 4 order. It is important to note the importance of having sufficient resolution (either through h or p) in the initial mesh in order to distinguish shocks, i.e. actual discontinuities, from smooth high-gradient regions when looking at high discontinuity sensor values.
We first run the solver on the initial mesh to obtain a base solution. We impose slip wall boundary conditions (BC) on the surface of the profile and far-field BC at the external boundaries of the domain. We use the HLLC Riemann solver as described by Toro (2009). For the artificial viscosity, we tuned the solver parameters to = −1.2, = 0.7 and 0 = 1.0. Fig. 2c (left) shows that large values of the sensor are obtained in both shock areas but also near the leading and trailing edges. However, Fig. 2d (left) shows that artificial viscosity only triggers in the vicinity of the two shocks, proving adequate tuning of the artificial viscosity parameters. As a result, the simulation is stable and we are able to reach steady state. The flow solution in Fig. 2b (left) displays very thick shocks as expected on this relatively coarse mesh. We can also observe oscillations in the field past the strong shock caused by the under-resolution of the shock and the generation of entropy.

r-adaptation
From the base solution, we follow the workflow explained in Sect. 6. We first extract the shock regions: these correspond to the elements of non-null artificial viscosity in Fig. 2d (left). To these regions, we assign a shrinking factor of 0.1 and run the r-adaptation procedure. We obtain a new mesh which, for quality considerations, we re-optimise before simulation. The new mesh shown in Fig. 2a (right) shows refinement in the shock areas and consequently a slight coarsening outside of those zones. Shrinking is also observed, to a smaller extent, in the vertical direction due to the isotropy of the r-adaptation approach. However, the resulting mesh is clearly anisotropic and aligned to the presence of the shock.
We now interpolate the old solution onto the adapted mesh and run the simulation again. In order to avoid any instability of the solver due to the interpolation of the under-resolved shock onto the new mesh, we first run the solver over a few hundred time steps with a decreased step size. We then run the simulation, using the exact same artificial viscosity parameters, until steadiness is achieved. The flow solution in Fig. 2b (right) shows better resolution of both shocks as seen by the sharpness of the shocks. We also observe reduced oscillations in the wake of the strong shock. Figs. 2c & 2d (right) finally show that discontinuity, as per the sensor, is now observed in a narrower area and that the artificial viscosity reaches lower values.
The improvement of the resolution of the shock can be better seen in 2D plots. Fig 3 shows the Mach number on the surface of the profile. Because elements are so large in the initial mesh, the solution shows strong oscillatory behaviour, known as the Gibbs phenomenon described in Sect. 1. After r-adaptation, elements in the region of the shock are much smaller and, although oscillations are still present, they have both a smaller amplitude and a narrower range. This confirms the qualitative observation of the increased sharpness of the shock seen in Fig. 2b. Because there is less mesh movement at the weak shock, the reduction in the Gibbs phenomenon is also smaller.
We also want to draw the reader's attention to the importance of CAD. To retain accurate boundary representation, it is important that the r-adaptation code has access to a CAD system. In this instance, Turner et al. (2016) implemented NekMesh to use the Open Cascade, developed by Open Cascade SAS (2019), kernel as its CAD engine under a small wrapper layer. This allows NekMesh to query the geometry and ensure that all nodes remain on the boundaries at all time. This is what we call CAD sliding and we show this capability in Fig. 4 where nodes remain on the aerofoil surface throughout the r-adaptation process. Fig. 4 also shows that the optimiser is able to move nodes across large distances, as seen through the row of coloured elements before (left) and after (right) r-adaptation.

p-adaptation
After better resolving the shocks, we can now apply local p-adaptation for the smooth field. For this test case, we compare three scenarios. In the first scenario, we apply local p-adaptation without any restriction (see Fig. 6a) while, in the other two scenarios, we restrict the local polynomial order inside the shock areas. In the second scenario, we preserve the local polynomial order of the uniform = 4 order simulation of Sect. 7.1.1 (see Fig. 6b). In the third and last scenario, we decrease the local polynomial order inside the shock areas to the lowest user-allowed order (see Fig. 6c).
For these tests, we start from the field obtained at = 4 in Sect. 7.1.1 and use values of min = 2 and max = 6. The sensor is based on the density field and the solver default values of lower and upper sensor tolerances are respectively 10 −8 and 10 −6 . Fig. 6 shows the results. The figures on the left show a final map of the local number of modes (= + 1) after a steady solution is reached and, by extension, when the local polynomial order remains constant throughout p-adaptation steps. The number of DOF for each simulation is shown in Table 1. All scenarios produce fewer DOF than the simulation at uniform p, thanks to local p-coarsening in low-error regions. As expected, the unrestricted p-adaptation scenario increases the local polynomial orders of elements in the shock thickness to the maximum user-allowed value. This leads to a higher global number of DOF than the other two scenarios. Then follows the second scenario while the last scenario has the smallest number of DOF. Each of these DOF counts also translates into similar increases or decreases in compute times.
We compare these solutions to a reference solution computed on a very fine mesh. To evaluate the performance of each mesh and p-adaptation scenario, we look at the Mach number distribution on the surface of the aerofoil. We use the 2 -norm of the error, defined as is the Mach number of the test solution, is the Mach number of the reference solution and is the chord. Results are reported in Table 1. We first note that r-adaptation alone provides an important boost in terms of accuracy. Scenarios #1 and #2 both suffer a loss of accuracy due to the coarsening of the solution in large parts of the domain. This slight increase in the error, however, allows us to cut the number of DOF in half. Scenario #3, on the other hand, performs very poorly, with the error going even higher than on the original mesh. Decreasing the polynomial order inside the shock -a rather small region -allows us to save a few more DOF but at too great a cost. Fig. 5 shows a comparison of the Mach number (left) and artificial viscosity (right) fields for the uniform simulation and the three test scenarios. We observe little difference between scenarios #1 and #2. Scenario #3 on the other hand exhibits underresolution of the shock, seen through its thicker profile. This is consistent with the local element size and the lack of DOF in the thickness of the shock at lower order. As a result, the last lower-order scenario exhibits some oscillations in the wake, due to the generated entropy in the shock area. We can also observe that lower-order scenarios induce more artificial viscosity. This phenomenon is consistent with the previous assessment of the lack of resolution of the shock. The discontinuity sensor detects a certain lack of resolution and therefore more artificial viscosity is added to the system.  Overall all scenarios expectedly exhibit similar distributions of local polynomial order in the smooth field regions in Fig. 6. When analysing the distribution of local polynomial orders, we observe higher orders in the area above the strong shock and below the weak shock, in all scenarios. These areas were not detected in Sect. 7.1.1 as part of the shock due to the then underresolved and therefore too short shocks. Now that the shocks are better resolved, they reach further out and require additional resolution, in the form of higher polynomial order in this case. We also observe, in the lower-order scenarios, that parasite higher-order zones are created. This is especially obvious around the weak shock in the third scenario. This is due, as we noted above, to the thicker shock profile and therefore the need to add resolution around the shock. Fig. 6c (left) is consistent with this explanation as we observe a larger area of high sensor values, extending beyond the shock areas determined in Sect. 7.1.1.

Supersonic intake
This section illustrates the new approach on a test case with a more complicated shock pattern. The test case is that of a supersonic intake at ∞ = 3.0 first studied experimentally by Anderson, Wong, and Field (1970) and later numerically by Jain and Mittal (2006). The intake consists of two straight ramps inclined with respect to the incoming free-stream flow at angles of 7 • and FIGURE 4 Magnified view of boundary elements with CAD sliding enabled nodes for the NACA 0012 test case: original (left) and adapted (right) meshes. Three rows of elements have been coloured to demonstrate the large movement of nodes. 14 • respectively. The first ramp creates an oblique shock which impinges on the horizontal cowl and in turn leads to a complex pattern of reflecting oblique shocks throughout the diffuser of the intake. The difficulty here is the presence of multiple shocks with different orientations in the very narrow regions of the diffuser. We discretise the domain uniformly in the streamwise direction. We set an element size of 0.01 ( being the length of the intake) inside the intake and let it coarsen outside the intake up to an element size of 0.05 in the far-field. The mesh is curvilinear of order = 4 and it is optimised in the throat. Fig. 7a (left) shows what the mesh looks like inside the intake and in its immediate surrounding.
We run the solver at uniform order = 3 on the initial mesh to obtain a base solution. We impose wall BC on the surfaces of the intake, at the intake outlet BC we set a low enough pressure until a fully supersonic field is obtained ( = 0.9 inf ) and far-field BC at the external boundaries of the domain. We use Roe's approximate Riemann solver as described by Toro (2009). For the artificial viscosity, we tuned the solver parameters to = 0.0, = 0.0 and 0 = 0.1. Fig. 7c (left) shows that large values of the sensor are obtained in all shocks and that moderate values are obtained everywhere after the first upstream shock. However, artificial viscosity only really triggers in the vicinity of the shocks as Fig. 7d (left) shows, proving adequate tuning of the artificial viscosity parameters. We are able to reach a steady state solution thanks to a stable simulation. Just like for the NACA 0012 test case, the shocks demonstrate a thick profile, as can be seen in Fig. 7b (left), due to the relatively coarse local mesh as well as some oscillations near the leading edge of the cowl.

rr-adaptation
We again follow the workflow laid out in Sect. 6 except that we decide to run two rounds of r-adaptation. Each round uses a less aggressive shrinking factor of 0.5. Before each simulation, we again optimise the mesh for high-order mesh quality reasons.
The new mesh after one round of r-adaptation is shown in Fig. 7a (right). We observe refinement in all areas of interest and note that refinement is stronger in the area of the first upstream shock. Indeed, elements in the first shock are able to pull DOF from the freestream areas whereas elements inside the intake are interacting with each other. Refinement is nonetheless obtained in all shock areas and anisotropy naturally appears such that elements are shrunk in mostly the shock normal direction. The radaptation strategy works by pulling nodes together. Because the shrinking areas are long and narrow, nodes are naturally moved normally rather than tangentially to the underlying shock, without the need for the optimiser to be aware of the shock structures. We now run the solver on the new adapted mesh using the same solver parameters. A stable flow solution is obtained and shown in Fig. 7b (right). All shocks now appear sharper and the oscillations observed near the leading edge of the cowl has disappeared. Figs. 7c & 7d (right) also show that discontinuity, as per the sensor, occurs in a narrower region.
We then apply a second round of r-adaptation in the exact similar fashion: we isolate shock areas and use them as input for the optimiser. Fig. 8a (right) shows the final adapted mesh which demonstrates further refinement of the shock regions. We also notice that the oblique shocks inside the intake past the throat have moved upstream due to the refinement of the oblique shocks located upstream of the throat. While the r-adapted mesh could not capture these downstream shocks, the rr-adapted mesh can. By using a two-step approach, we are also able to pull more mesh nodes together than when using a one-step approach.
This becomes even more obvious when looking at the plot of the Mach number in Fig. 9. The initial mesh largely overestimates the values of the Mach number past the throat ( ∕ ≈ 0.57). This is mostly solved by r-adaptation although further improvement is obtained through rr-adaptation-adaptation. This is due to the good resolution of upstream shocks through the clustering of DOF.

p-adaptation
We now apply p-adaptation to the rr-adapted mesh. For this test case, we focus on unrestricted p-adaptation where the local polynomial order inside elements is left free to change, even in shock areas. We start from the field obtained at = 4 in Sect. 7.2.1 and use values of min = 2 and max = 6. We again use a sensor based on the density and solver default values for the thresholds.
First, we observe that no steady state is achieved. Upon inspection, we notice that the system jumps back and forth between two states at each p-adaptation cycle. The two states correspond roughly to coarser and finer resolved fields. In the coarser resolved state, sensor values in shock areas are high. At the end of the p-adaptation cycle, these large sensor values trigger an increase in local polynomial order of a number of elements. Simulation goes on and the finer resolved state is obtained where sensor values are low. This in turn triggers a decrease in local polynomial order of the same elements, returning the system to the former coarser resolved field. This is shown in Fig. 10 with the coarser resolved state on the left and the finer one on the right. We explain this behaviour by a naive p-adaptation approach using simple sensor thresholds. The problem is highly nonlinear and non-local and error from refining/coarsening regions propagates along characteristics. The non-adjoint nature of the refinement strategy is bound to produce this sort of behaviour.
Nevertheless we observe that additional resolution in the form of higher local polynomials is found in sensible areas: in the shocks, inside the intake (especially in the throat) and right above the coil. The only very high polynomial orders are obtained in the shocks whereas smooth regions reach order = 3 at most. The number of DOF for each simulation and state is shown in Table 2. Referring back to Fig.9, we can see that little difference in the solution appears from rr-adaptation to rrp-adaptation despite the decrease in number of DOF.

CONCLUSIONS
In this paper, we have presented a novel strategy for adaptive simulations, based on a combination of both -and p-adaptation. The proof-of-concept work applied here takes advantage of both strategies in different manners, as appropriate for the simulation of compressible flows containing shocks. We achieve mesh movement required for r-adaptation through the use of a variational optimisation strategy, using the combination of a local discontinuity sensor and a target element size in order to effectively cluster DOF in the presence of shocks and more sharply simulate their features. At the same time, we apply a p-adaptation technique in the rest of the domain in order to benefit from the spectral rate of convergence of high-order discretisations for smooth solutions. The simulation is effectively stabilised through the use of an artifical diffusion term, again using the local discontinuity sensor.        The proposed strategy exhibits a number of benefits from a computational perspective, as seen in the results presented in the previous section, where the canoncial NACA 0012 test case and a more challenging supersonic intake have been examined. The main benefit of this dual-adaptive technique is that we are able to significantly reduce the number of DOF required to resolve a given simulation, when compared to a uniformly refined grid or using solely r-adaptation. Table 1 shows that, for the various -refinement strategies considered, the error when compared to a very fine solution remains roughly the same, whilst the simulation requires only 50% of the DOF of the original simulation. This has important consequences from the perspective of computational efficiency, since a significant reduction in the number of DOF will lead to a reduction in runtimes. Likewise, the cost of operations per DOF is reduced as the polynomial order decreases, which offers the opportunity to further reduce computational cost. The rp-adaptation technique therefore permits an effective balance to be achieved between the attained error and simulation expense.
From the context of more general conclusions of our results, we demonstrate that care must be taken when selecting a padaptation strategy. In particular, the NACA 0012 simulations demonstrate that p-coarsening can have important negative effects on the solution for minimal computational gains. Additionally, the supersonic intake exhibits a complex shock pattern. Because of the complexity and strength of the reflecting shocks, we show that multiple r-adaptation steps are not only possible but

FIGURE 10
Comparison of the fields for the intake in its coarser (left) and finer (right) resolved states during unrestricted p-adaptation.
desirable. Despite the lack of nodes to redistribute inside the intake, sufficient mesh deformation is achieved to better capture the different shocks.
Although the overall strategy has been shown to be effective, it is important to emphasise that some of the benefits we highlight in this work can be attributed to our particular implementation of the r-adaptation technique. In particular, the use of the variational framework yields several advantages. Firstly, the use of a target element size allows the grid to deform in an anisotropic manner within restricted regions of the domain. Even when the deformation is substantial, this still permits a valid grid to be obtained, as shown in Fig. 2. Secondly, the ability to conform to complex CAD surfaces and curves whilst permitting nodes to slide across them is clearly important in the context of this work, where shocks arise at or near solid surfaces. This functionality can be difficult to achieve in other mesh deformation techniques, particularly those that require the solution of a PDE system of an appropriate solid body model.
Finally, we note that there are a number of clear directions for potential future work in this area. An extension of this method to transient flows, especially with moving shocks, would constitute an interesting application of this rp-adaptation strategy. The variational moving mesh framework would be able to track shocks throughout the simulation without the need to generate a new mesh. With preserved mesh connectivities, the system of equations would not need to be re-built at each adaptation step. This is especially desirable on large meshes and large HPC-based simulations where I/O and inter-node communication can incur significant expense.