A free–energy stable p–adaptive nodal discontinuous Galerkin for the Cahn–Hilliard equation

A novel free–energy stable discontinuous Galerkin method is developed for the Cahn– Hilliard equation with non–conforming elements. This work focuses on dynamic polynomial adaption (p–reﬁnement) and constitutes an extension of the method developed by Manzanero et al. in Journal of Computational Physics 403:109072, 2020, which makes use of the summation–by–parts simultaneous–approximation term technique along with Gauss–Lobatto points and the Bassi–Rebay 1 (BR1) scheme. The BR1 numerical ﬂux accommodates non– conforming elements, which are connected through the mortar method. The scheme has been analytically proven to retain its free–energy stability when transitioning to non-conforming elements. Furthermore, a methodology to perform the adaption is introduced based on the knowledge of the location of the interface between phases. The adaption methodology is tested for its accuracy and eﬀectiveness through a series of steady and unsteady test cases. We solve a steady one–dimensional interface test case to initially examine the accuracy of the adaption. Furthermore, we study the formation of a static bubble in two dimensions and verify that the accuracy of the solver is maintained while the degrees of freedom decrease to less than half compared to the uniform solution. Lastly, we examine an unsteady case such as the spinodal decomposition and show that the same results for the free–energy are recovered with a 35% reduction of the degrees of freedom for the two–dimensional case considered and a 48% reduction for the three–dimensional case.


Introduction
A variety of different approaches and methodologies have been developed throughout the years to model and simulate systems of two or more immiscible fluids. These methods can be divided into two large categories. The interface-tracking and the interface-capturing techniques. The former includes the Marker-And-Cell (MAC) and front-tracking methods whereas the latter include the widely known Volume Of Fluid (VOF), level-set and phase field methods. A review of the characteristics and features of the various methods can be found in [1].
In this work, we focus on phase-field methods and especially the Cahn-Hilliard model [2]. Its use for multiphase modeling has growing popularity in the recent years because of the favorable total mass conservation characteristics [3] as well as the bounded behavior of the free-energy [4].
The essence of using a diffuse interface method such as the Cahn-Hilliard model is that the sharp interface between the two immiscible fluids is represented by an interface of a finite width throughout which, the fluid thermodynamic properties smoothly vary from those of one phase to another. Typically, the physical width of the interface is of the order of nanometers [5]. However, this constitutes an unrealistic target as to numerically capture this interface would require a large amount of resources even for the simplest of the cases [6]. The answer to this is partially derived from the sharp interface limit of the Cahn-Hilliard model. Through asymptotic analysis [7], it can be proved that the Cahn-Hilliard model converges to the physical solution even with the use of a substantially larger interface. This characteristic makes a simulation cost-effective while retaining the desired accuracy.
The discretization scheme adopted is the nodal Discontinuous Galerkin Spectral Element Method (DGSEM) based on previous work presented by Manzanero et al. [8]. The DGSEM offers great flexibility as the solution is represented through an arbitrary approximation polynomial order and supports the use of unstructured meshes of curvilinear hexahedral elements to approximate complex geometries. Furthermore, another advantage of the DGSEM exploited in this work is that the approximation order can vary across elements. There is a limited amount of publications concerning the use of discontinuous Galerkin methods in order to solve the Cahn-Hilliard equation [8][9][10][11][12][13][14] and there is still a variety of aspects regarding the efficiency, robustness and accuracy that should be addressed.
The DGSEM version of this work uses Gauss-Lobatto (GL) points, which through the Summation-By-Parts Simultaneous-Approximation Term (SBP-SAT) property [15,16] allow the derivation of free-energy stable schemes. There is a plethora of work focused on the construction of entropy and energy stable schemes for the discontinuous Galerkin method [15,[17][18][19][20][21][22][23] and the references therein. A review of entropy stable DGSEM schemes is given in [24]. As for the Cahn-Hilliard equation, a DGSEM with the SBP-SAT property has been developed in [8,9]. These schemes however, are designed for conforming elements and they should be modified accordingly such that they can retain their stability properties when using non-conforming elements.
There is a growing popularity of entropy stable schemes with special consideration on the extension of their characteristics to non-conforming elements. These ideas can contribute to the creation of a free-energy stable scheme as the one developed in this work. The foundation for the establishment of this family of methods has been laid in [25][26][27][28]. One of the first works is the one presented by Carpenter et al. [29] which makes a first attempt to create an entropy-stable scheme using the local discontinuous Galerkin method. Again, the scheme uses the SBP-SAT property of the Gauss-Lobatto (GL) points and is proved to be entropy-stable for the compressible Euler equations for unstructured hexahedral meshes. This method has been extended to the Navier-Stokes equations by Parsani et al. [30,31]. Then in [32] the classic mortar method is modified to guarantee entropy stability of the discontinuous Galerkin method. Following that work, the authors in [33] extended the scheme from the compressible Euler equations to the compressible Navier-Stokes equations and into the use of curvilinear elements. In addition, [34] is also an extension of the aforementioned work and focuses on an entropy stable discontinuous Galerkin scheme. Another work on non-conforming curvilinear elements for hyperbolic conservation laws is presented in [35]. All the literature presented so far is based on the Euler and Navier-Stokes equations and up to our knowledge, an equivalent free-energy stable scheme for non-conforming elements has not been previously established for the Cahn-Hilliard model.
For diffuse interface methods, such as the Cahn-Hilliard, the region of interest is that of the interface. More specifically, within the interface there are steep gradients, which should be adequately resolved. On the contrary, on the bulk of each phase, the phase-field variable is constant and the resolution can be lowered. Thus, the resolution can be coarsened in the regions away from the interface and refined within the vicinity of the interface. The interested reader can find a wide array of publications that address the aforementioned issue in the context of mesh adaption (h-refinement) [36][37][38][39][40][41][42][43][44] with respect to the refinement indicator, mesh refinement strategy and dynamic adaption. These correspond to different frameworks, employing numerical schemes such as finite differences [38,39], finite volume [40,41] and finite elements [36,37,43,44]. Also, there are some works on finite elements that identify the effect of the polynomial order on the accuracy of the solution [13,45], but polynomial adaption (p-refinement) is not addressed. Up to the authors' knowledge, there is no similar work performed on the effectiveness of local polynomial adaption for phase-field problems when used in conjunction with the discontinuous Galerkin scheme.
There are various methods that can be used to mark the refinement region for Adaptive Mesh Refinement (AMR), which could potentially be applied to polynomial refinement. Similar studies on dynamic mesh adaption for diffuse interface models have been performed with the use of gradient based indicators [37,43], Legendre discretisation spectrum extrapolation [46], summation of the tail of the spectrum [47], flux jump of the Laplace operator [36], the location of the interface [38,48], as well as other error indicators [49][50][51]. Additionally, there is a theoretical study on the a priori determination of the error for the discontinuous Galerkin method when used for the Cahn-Hilliard equation [52]. A possible alternative approach for the polynomial refinement of the Cahn-Hilliard is the truncation error estimation as presented in [53] and applied to the Navier-Stokes equations in [54][55][56], which could allow anisotropic refinement in each element. In this work, we take advantage of the knowledge of the location of the interface to perform the adaption.
The aim of this work is to create a scheme that retains the accuracy and stability characteristics when transitioning to the use of the non-uniform polynomial order across the domain. To achieve this, the first step was to modify the BR1 numerical flux to handle such elements and incorporate the mortar method [57]. The modified scheme is then proved analytically to be free-energy stable for the general case of non-conforming element boundaries following the methodology in [8]. The time marching is performed through a first order IMplicit-EXplicit (IMEX) scheme. This choice has been made to alleviate the stiffness of the system due to the higher order derivatives. The free-energy of the system is proved to be bounded for the continuous and the discrete time settings.
We also develop a methodology to perform automatic p-adaption for the Cahn-Hilliard equation when using the DGSEM. This is based on the location of the interface and it takes advantage of the attributes of the scheme's spatial locality. The DGSEM offers the opportunity to augment the polynomial order within each element of the mesh and also individually in each direction. Therefore, only the elements that contain part of the interface are refined while the rest are coarsened. Thus, through this simple and efficient process we can achieve a reduction in the degrees of freedom without further processing of the mesh. The accuracy and convergence characteristics are subsequently examined through a series of different tests and direct comparison with the conforming version of the scheme.
The rest of the work is organized as follows. First the model used and its characteristics are specified in Sec. 2. Then the spatial and temporal discretisation methods are introduced in Sec. 3, as well as the modified BR1 numerical fluxes that allow non-uniform polynomial order in Sec. 3.2. The adaption methodology is described in Sec. 4. In Sec. 5 the foundation of the method is laid through an analytical proof of the free-energy boundedness of the scheme. In Sec. 6 three different test cases are presented along with the results from the developed method and a comparison is made with the conforming solver version with respect to the achieved accuracy and degrees of freedom.

Cahn-Hilliard equation and continuous energy estimates
This work is based on the Cahn-Hilliard equation, which can be utilized to describe a vast variety of physical phenomena such as the dynamics of phase separation of two and N phase flows, binary alloys, tumor growth, etc. [3]. The application of interest in this particular case is the study of two phase flows. The model consists of a constant mobility parameter and a polynomic double-well chemical free-energy [58,59]. The Cahn-Hilliard equation is defined as The phase field variable φ represents the concentration of each phase and satisfies (1), M represents the mobility which is a positive parameter, Ω is the physical domain (with boundary ∂Ω), and w is a scalar field representing the chemical potential. The chemical potential is designed to minimize an arbitrary free-energy functional, F(φ, ∇φ), which depends on the phase field and its gradient, A homogeneous Neumann boundary condition is applied for the chemical potential in order to ensure that mass is conserved [60], The free-energy consists of two terms that impose opposing effects. The chemical free-energy, ψ, which drives phase separation, and the interfacial energy 1 2 k|∇φ| 2 , which promotes homogenization by penalizing the existence of gradients, which consequently translates into the existence of an interface, In (4), F v (φ) and F s (φ) represent the volumetric and surface free-energies accordingly. The term g(φ) represents a boundary energy that will also be minimized with appropriate boundary conditions, and k is the interfacial energy coefficient. The minimization is performed by linearization of the free-energy (4) around an equilibrium solution, where δφ acts as a small perturbation. Since we will also apply Neumann boundary conditions for φ, the perturbation δφ is not restricted to vanish at the boundaries ∂Ω. The integration of the second term of the first integral in (5) by parts, yields both the chemical potential definition, and the appropriate Neumann boundary conditions prescription, In this work we use a polynomial double-well function for the chemical free-energy [2], and a linear function for the boundary free-energy, since they represent standard choices in the literature, but other choices that are not covered here exist (e.g. logarithmic chemical free-energy [61]). An exact solution of the one-dimensional steady state Cahn-Hilliard equation (1), with the chemical free-energy (9) and the free-energy (4), in an infinite domain with φ (±∞) = ±1 is given by To facilitate the stability analysis of the system as well as the discretization, and following the processes highlighted in [8,10], the Cahn-Hilliard equation (1) is transformed into a system of four first order equations, To create the weak form, we have introduced two auxiliary variables q = ∇φ and f = ∇w. We multiply (12a) and (12c) with the arbitrary scalar test functions ϕ Φ and ϕ W respectively and (12b) and (12d) with the arbitrary vectorial counterparts ϕ F and ϕ Q . Then we integrate over the domain Ω and derive four weak forms, where the operator f, g is the L 2 inner product

Continuous free-energy stability
The analytical proof of the free-energy boundedness with a given initial condition is given in [8].
The process makes use of the transformed system (13) and then through some manipulation of the equations the following inequality can be derived, As a result, the Cahn-Hilliard equation (1) with chemical potential (7) and Neumann boundary conditions (3) and (8) guarantees that the free-energy F evaluated at any time instant T , as defined in (4), is bounded.This is the property to be mimicked by the subsequent approximation.

The nodal discontinuous Galerkin spectral element method
In this section we briefly introduce the underlying theory and the construction of the nonconforming nodal DGSEM. The interpolating nodes are the GL points as they satisfy the SBP-SAT property (19). This is a crucial aspect to prove the scheme's stability, avoiding the use of exact integration. Furthermore, a non-conforming DGSEM is constructed, where the polynomial order can vary across different elements [62].
The computational domain Ω is tessellated into non-overlapping hexahedral elements, which are then geometrically transformed from a reference element e = [−1, 1] 3 by means of a transfinite mapping. This mapping relates the physical ( x = (x 1 , x 2 , x 3 ) = (x, y, z)) and the local ( ξ = (ξ 1 , ξ 2 , ξ 3 ) = (ξ, η, ζ)) coordinates, The solutions and functions are approximated by order N e polynomials in an element e, where the approximation order N e can vary from element to element. In (17), l j are the Lagrange interpolating polynomials whose nodes are a set of Gauss-Lobatto points in the reference element e and U ijk (t) are the (time dependent) nodal values of an arbitrary function u. The notation is as follows: we use lower cases for the exact functions, whereas upper cases represent their polynomial approximation. We approximate the integrals with quadrature rules that use the same GL nodes as those that represent the solution, where w i are the quadrature weights [63]. This provides an exact integration for 2N e − 1 order polynomial (see [64]). The choice of the GL nodes is essential since it makes the quadrature rule (18) to satisfy the discrete SBP-SAT property, that is, a discrete Gauss law, which in one dimension is The SBP-SAT property makes it possible to resemble the continuous analysis that proves the boundedness of the free energy F in (15) discretely.
The transformations from the reference element to the physical space must create a watertight mesh (i.e. without gaps across the elements). This constraints the mapping functions at the inter-element faces to be the same, From the mapping, we compute a set of covariant and contravariant vector bases, The relation between the covariant and (volume weighted) contravariant bases is where J is the Jacobian of the transformation, J = a 1 · ( a 2 × a 3 ). The volume weighted contravariant bases satisfy the continuous metric identities, In the discrete setting, we must ensure that the approximation of the metrics is free-stream preserving. Having a watertight mesh is a necessary condition but not sufficient. For free-stream preservation, the mesh has to satisfy two additional conditions [65], 1. Condition (F): the projection of the discrete volume weighted contravariant bases at the faces has to be also continuous J a i e L f ace = J a i e R f ace . The symbol J represents the polynomial approximation of the Jacobian. The mapping function is approximated with the interpolation operator (17). However, although the mapping is represented by order N e polynomials, the genuine order of the mapping is given by the approximation order of the faces. We highlight the construction of a watertight mesh that satisfies the two conditions: • Edges: the approximation order of the edges in the mesh has to be unique, and the face functions that share an edge must reduce to the same curvilinear function at the edge. The order of an edge shared by various faces has to be N edge = min (N f 1 , N f 2 , ...) /2 for general three-dimensional non-conforming elements, and N edge = min (N f 1 , N f 2 , ...) for two-dimensional, two-dimensional extruded, and conforming problems, at most.
• Faces: the approximation order of the faces in the mesh has to be unique and it must be N f = min (N e L , N e R ) /2 for general three-dimensional non-conforming elements, and N f = min (N e L , N e R ) for two-dimensional, two-dimensional extruded, and conforming problems, at most.
• Volume: the contravariant basis have to be computed in a curl form [66] J We use the contravariant basis to transform the differential operators from physical to computational space. The divergence of a vector is where, is the contravariant flux, and The gradient of a scalar is transformed as [66] ∇u Lastly we approximate the three dimensional integrals in an element by GL tensor product quadratures which allows us to write the discrete Gauss law as in [67] ∇ ξ U,F In (30),n is the unit outward normal vector at the reference element faces, dŜ is the surface local integration variables (dŜ i = dξ j dξ k for i-oriented faces). To compute the surface integral, the two dimensional quadratures in each of the six faces that define the element are defined Moreover, surface integrals can be written in either physical or computational space. The relation between the two spaces is where we defined the face Jacobian J i f = J a i . An equivalent relation can be deducted for the surface flux in the reference element,F ·n, as well as the physical, F · n, through the relation Therefore, quadratures can be represented both in physical and computational spaces, and the use of one form over the other depends on whether we study an isolated element (computational space) or the whole combination of elements in the mesh (physical space).

The mortar element method
Following [57], we use two projection operators P lh and P hl from the order N l (low order) space to its N h > N l (high order) counterpart and viceversa known as the mortar method [46]. Let F ∈ P N l and G ∈ P N h two polynomial functions. Then, the projection operator is designed to satisfy A consequence of (35) is that the operator to augment the polynomial space differs to the opposite one, known as restriction. To augment the polynomial space one simply uses a interpolation operation, that is, we evaluate the Lagrange interpolating polynomials of the low-order N l space at the GL nodes of the high-order N h space. To fulfill the condition (35), the backward projection (restriction) from N h to N l must be where w Ne j are the quadrature weights for an element e with polynomial order N e . Constructed in such way, the operators are not invertible, that is P hl P lh = I h and P lh P hl = I l . This means that when an arbitrary polynomial A ∈ P N l from the lower space is projected to the mortar, it is not recovered by using the restriction operator, An exception to (38) is any constant function, P hl (P lh (k)) = k. Therefore, we can show that although the polynomials in (38) differ, the mortar is still able to keep the integral value,

Discontinuous Galerkin spectral element approximation of the Cahn-Hilliard equation
We now assemble the discrete version of (13). The first step is to transform (13) into the local coordinate system as described in (25) and (28), and to get the weak form of the system in the reference element E. To do that, we restrict the test functions ϕ Φ , ϕ W (scalar), ϕ F , ϕ Q (vectorial), to the order N e polynomial space, The following step is to integrate the right hand side terms that contain a ∇ ξ operator by parts, to replace the continuous functions with their polynomial approximations and to replace exact integrals by quadratures. Furthermore, we apply the discrete Gauss law (30) to (41b) and (41d).
The terms with star superscript in (41) are the numerical fluxes, which make the flux uniquely defined at the boundaries. In this work, we use the Bassi-Rebay 1 scheme (BR1) [68], which handles non-conforming interfaces. Without loss of generality, we consider the inter-element face with orders N l < N h . The mortar method computes the scalar fluxes on the higher order element, which are then transferred to the lower order element through the restriction operator, Whereas for the vector fluxes, we transform the contravariant fluxes from the lower order to the higher order element and then we take their difference (since the normal vectors are opposite) and add an interface stabilizing term, In (43), |J a| h n h and |J a| l n l are the (scaled) normal vectors at the face, and they satisfy |J a| h n h = − |J a| l n l = −P lh |J a| l n l (Condition (F)).
For Neumann boundary conditions, we use the adjacent element interior value to compute the gradients in (41b) and (41d), and we directly impose Neumann boundary values for divergence weak forms (41a) and (41c), Note that we have presented the physical interface fluxes, rather than the contravariant fluxes, so that we can relate in a straightforward manner the interface values shared by two elements (recall that the same physical flux yields different contravariant flux values, as it depends on each element geometry).

Temporal Discretization
The Cahn-Hilliard equation incorporates higher order derivatives which make the system stiff.
To combat the stiffness, a first order IMplicit-EXplicit (IMEX) scheme has been employed. Throughout this work a constant timestep ∆t has been used. Following the theory developed in [8,69], the scheme takes the following form: The notation follows that n will be used as a superscript to denote state values at t n = n∆t.
The chemical free-energy ψ (φ) is non-linear and will be treated explicitly as denoted in (46). The treatment of the interface energy term k∇ 2 φ is controlled by the K 0 parameter. For K 0 = 0 there is a fully explicit time marching, for K 0 = 1/2 we get the Crank-Nicolson scheme, whereas for K 0 = 1 it is a fully implicit scheme. For the numerical tests of Sec. 6, the fully implicit treatment has been utilized. The second term of the right hand side of (46) is an additional numerical stabilization according to [69]. It has been shown in [8] that a value of S 0 = 1 the system is stable for Φ ∈ [−1, 1] and thus has been chosen for the following numerical tests.
Applying the scheme (46) to the discretized scheme of (41), the fully discrete discontinuous Galerkin approximation is obtained The superscript θ has been used for variables (e.g. F θ or W θ ) that are not directly evaluated at t n or t n+1 with the IMEX strategy, but on a combination of those depending on the different terms involved in (47c).

Heuristic p-adaption methodology
In this section we describe a methodology to automatically adapt the polynomial order for solutions of the Cahn-Hilliard equation. This methodology exploits the characteristics of the introduced numerical method which permit the use of different order of accuracy in each element. The indicator identifies the elements that contain at least one point within the interface region. The latter has been defined as the region where the phase field parameter φ ranges from −0.9 < φ < 0.9 and is a choice which is usual in diffuse interface literature [70]. If an element is marked for refinement, then it is refined to a user specified polynomial order in all directions. The adaption process is presented in Algorithm 1.
To enhance the robustness of the method, and since the movement of the interface cannot be predicted when solving the Cahn-Hilliard equation, a buffer region of refined elements is applied to all the neighboring elements of those that contain part of the interface, as presented in Algorithm 1. On the contrary, the elements that are not marked, are coarsened also to a user specified level. In this work we have used a coarse level of N Coarse = 2. Simulations have also been conducted with a coarse polynomial order of N Coarse = 1 and we noticed that spurious oscillations were introduced in the free-energy and thus this option has been omitted.
In order to avoid steep changes in the polynomial order across the elements, two different jump restrictions have been applied to bridge the fine and coarse element levels. In both cases, the polynomial order of the elements in fine level remains intact, and modifications take place in the neighboring coarse level elements. The first criterion restricts the polynomial order jump between two subsequent elements not to be greater than unity and from this point on, we will refer to that as the N − 1 condition. The other criterion is less restrictive and dictates that for two adjacent elements with N l < N h , the N l will take a value of N l = 2N h /3. The latter will be referred to as the 2N/3 criterion for the rest of this text. These conditions have been derived from the work of [62].
One of the most important benefits of this adaption approach, is that the computational cost is kept to a minimum. The method, calculates the new polynomial order for each element and performs the interpolation or projection from the previous to the following polynomial orders. In this case we use the operator P P N (35) to perform an L 2 projection from order P to order N as presented in [63], This method is effective since the solutions of the Cahn-Hilliard equation typically involve a single scale, that of the interface. This flexibility also stems from the use of a higher order method such as the discontinuous Galerkin. Due to the higher order approximation, the element size can be substantially larger compared to classic lower order methods and thus the issues for the interface-marking approach detailed in [50] such as the refinement or coarsening of the area around the interface are bypassed. In addition, for phase field simulations, the target is to have a sufficient number of solution points within the interface in order to capture it appropriately. For other discretisation methods, this number is defined to be at least 5 points within this region [71], while some researchers prefer a number close to 10 [72,73]. In this work, we use a varying number of points within the interface that ranges from 3 to 30 to assess the convergence characteristics of the scheme.
Therefore, since we have a-priori knowledge of the element size, the interface width and the desired value for points within the interface, it is a very straightforward and low-cost method which does not require a rigorous tuning of parameters such as in various marker methods used for AMR [37,43,50]. Thus, as presented in Algorithm 1, the user specifies a value for the maximum and the minimum polynomial order to be used, the iteration interval upon which the adaption will take place and lastly the criterion for the polynomial order jump. if mod iter f adpt = 0 or iter = 1 then 8: for i ← 1,number_of_elements do 9: if any |Element(i).φ (x, y, z)) | ≤ 0.9 then

Stability analysis
In this work we focus on the p-non-conforming extension of the stability analysis methodology presented in [8]. The initial steps of the analysis are briefly described and a more detailed description can be found in [8]. This analysis focuses specifically on the treatment of the interelement boundary terms for elements with different polynomial orders to prove that the scheme is free-energy stable.

The time derivative of the free-energy is identified from the terms
This results in the following equation for an element e, where IBT and PBT are the interior and physical boundary terms. As derived in [8], the PBT represent the surface free-energy, thus, the equation for the total free-energyF = e F E,Ne + boundary faces

∂E,Ne
GdS is A stable approximation has IBT 0. It is in the interior boundary terms where the stability of the inter-element coupling through the mortar element method is assessed. Each interior face has the contribution of the two adjacent elements with orders N l < N h , For simplicity, we first study the second and fourth terms of (53) (the part that involves W and F ). We replace the inter-element fluxes (42) and (43), and we apply property (35) to the second term in (53), We add the (negative) contributions from the fourth term in (53) and the transformed lower order integral of (55), We proceed similarly with the first and third integrals of (53), for which we use the same numerical fluxes, which gives, Finally, following [8] we define an augmented free-energyF σ that includes the interface penalization, and that satisfies the free-energy equation, The physical dissipation and the numerical dissipation are both responsible for decreases in the free-energy. The numerical interface penalization is not needed for stability, hence a valid scheme can have σ = 0, but having σ > 0 usually enhances the accuracy of the solutions [8,74].

Stability analysis of the non-conforming fully discrete system using the IMEX integrator
To begin the stability analysis for the fully discrete system the following steps are carried out in accordance to [8]. The analytical relations can be found therein, as we will mainly focus on the effect of having p-non-conforming elements which are connected through the mortar method as described in Section 3.1. The initial steps for the stability analysis are briefly provided: 1. We evaluate (47d) at t n+1 and t n and subtract them. Then we divide the by ∆t and set φ Q = Q n+1 .
4. We then subtract the equations derived in Steps 2 and 3 and multiply by ∆t .
5. The following step is to express the chemical free potential ψ (φ) through a Taylor expansion using the definition from (9).
6. Through some manipulation, which is described in detail in [8], we end up to an expression for the discrete volumetric free-energy F n,E,N v within each element, which is defined as F n,E,N v = J Ψ n + 1 2 k Q n · Q n , 1 .
7. The last step is to sum over all the elements within the domain Ω.
Following the aforementioned steps we derive (61) for the discrete volumetric free energy: where IBT represents the inter-element coupling terms of non-conforming elements, with approximation orders N l and N h respectively, as defined in (62). The terms PBT arise from the exterior boundaries (66) and the application of boundary conditions. The last term is the numerical dissipation of the IMEX scheme, which is negative, and there is a detailed proof in [8] about its contribution in the decrease of the free energy. As the non-conformity of the domain does not affect or alter that term, it will not be presented analytically. Lastly, the first term of the right hand side of (61) will always be negative and represents the physical dissipation arising from the interior of each element. Following the aforementioned steps, the interior boundary terms can be expressed as in (62). This is the equivalent version of (53) from the continuous stability analysis: where ∆Φ = Φ n+1 − Φ n . The second and fourth terms from (62), and following the same procedure as in (56), are equal to: The first and third terms of (62), following the process of (57), are equal to: We manipulate (64) by adding and subtracting Φ n and we arrange it to get The boundary terms remain unchanged because we use same approximation in the interior as the exterior for every boundary: Gathering the results of (63), (65) and (66) and substituting in (61), we derive (67) for the volumetric free-energy. We define an augmented free-energy F n,N,σ Equation (67) confirms that the free-energy remains bounded in the discrete time setting, as the four terms on the right hand side of (67) are dissipative.

Numerical results
In this section, we consider three different numerical experiments: a one-dimensional interface test case, the formation of a circular static bubble and that of a spinodal decomposition in two and three dimensions. The one-dimensional test case is used to initially evaluate the accuracy of the adaption compared to the conforming version of the scheme as well as to quantify the error introduced when adapting along the direction of the interface. The static bubble formation is an unsteady case used to assess the accuracy of the scheme in two-dimensions and showcase the ability to achieve similar level of accuracy through a significant reduction in the degrees of freedom. Lastly, the spinodal decomposition is a classic problem of the Cahn-Hilliard equation.
We examine a two-and a three-dimensional test case of a spinodal decomposition. We show that the same results for the free-energy and final state can be recovered with the use of adaption and present the course of the reduction of the degrees of freedom as the solution evolves. These test cases have been chosen to illustrate the capabilities, robustness and accuracy of the scheme.

One-dimensional interface
The first test case to be examined is an one-dimensional interface between two immiscible fluids.
It is a case that proves the effectiveness of this methodology and can be used as a preliminary study to assess the characteristics of the scheme. Moreover, the impact of the polynomial order jump condition on the accuracy among neighboring elements is assessed. The initial condition is a sharp jump of the phase-field variable φ between the two phases is defined as The Neumann boundary condition (8) is enforced at both ends of the domain. The model parameters used in this case are presented in Table 1. Recall that the parameters K 0 and S 0 are related to the IMEX time marching scheme. The backward difference Euler has been chosen (K 0 = 1) as well as the value for S 0 = 1 that guarantees non-linear stability [8].   Table 2: mesh specification for the test case of the steady state one-dimensional interface and relative element size compared to the interface width as approximated by (11).
As the simulation evolves, the system converges to a steady state. The exact solution of the steady Cahn-Hilliard equation (1) with the chemical free-energy (9) is (11). We chose to compare the convergence of the scheme based on a very fine solution of N = 10 and not with the analytical solution due to the domain not being extended to infinity. The initial condition of the problem along with the final solution and the analytical solution of (11) are presented in Figure 1. The adaption interval in this test case does not correlate with the accuracy of the results, as this is a steady case and the adaption performed on the first timestep is sufficient. We first establish the convergence rates and characteristics when uniformly refining the polynomial order of the solution, as presented in Figure 2. More specifically, the convergence test has been performed for three different meshes. The coarser case corresponds to a cell size which is equal to the interface width, as approximated from the steady state solution. The subsequent two meshes have cell sizes of half and a quarter of the initial cell size respectively. From the results presented in Figure 2, we confirm that the polynomial refinement leads to exponential convergence rates. The convergence rate has an even-odd behavior which is also present in the manufactured solution case in [8]. Through this test case, the effect of the polynomial refinement has been showcased as well as the dependency on the cell size of the mesh used.  Table 2.
Having established for this case that the polynomial refinement leads to exponential convergence, we validate that these characteristics are preserved when performing polynomial local refinement and coarsening. In all cases, the initial setup has been a uniform fine level polynomial order across the domain. This polynomial order is retained in the region of the interface, whereas there is coarsening in the bulk of the phases. The lowest polynomial order chosen for the bulk of the phases is N Coarse = 2 as described in Sec. 4. This case also incorporates an initial test for the jump of the polynomial order across different elements. Since it is a one-dimensional case it does not unveil the effect and the error that might occur due to the use of the mortar method. However, it has been used as a mean of identifying errors that occur while changing the polynomial order in the direction of the interface. Two cases for the jump of the polynomial order between two subsequent elements have been tested, the N − 1 and the 2N/3 criteria as introduced in Sec. 4.
The results of the adaption process are presented in Figure 2. The results with the use of adaption match those of the uniform polynomial order for all meshes and polynomial orders tested. Furthermore, the reduction achieved in degrees of freedom throughout the simulation of the one-dimensional interface is presented in Table 3.

Circular static bubble
The second case to be studied is that of the formation of a static bubble. The initial condition is a square at the center of the domain with the phase-field parameter taking the value of φ = −1 and the rest of the domain the value of φ = 1, The initial condition is presented in Figure 3(a). Subsequently, the solution evolves to a state of minimum energy, which is a circular bubble as presented in Figure 3(b). The parameters of the simulation are presented in Table 4. The final solution time is t = 100. Two different Cartesian meshes have been considered. The first has an element size which is half of the interface width and the second being a quarter of the interface width. Also in this case, the element size was based on the knowledge of the parameter in the Cahn-Hilliard equation that controls the interface width. The specification of the mesh is presented in Table 5. All boundaries are treated with periodic boundary conditions. A suitable adaption interval for this test case and the meshes specified in Table 5 has been identified to be every 500 iterations. It has to be mentioned that this value is dependent upon the element size, the dynamics of the problem as well as the timestep chosen for each particular simulation.   Table 5: mesh specification for the test case of the static two-dimensional circular bubble and relative element size compared to the interface width as approximated by (11).
The results using the adaption algorithm on this problem compared to the uniform solution are presented in Figures 5 and 7. The data points for the adapted results in Figure 5 correspond to Mesh 1 and a fine level polynomial order of N Fine ∈ {2, 3, 4, 5, 6, 8} whereas N Coarse = 2. For this particular mesh, the results amongst the cases with the different polynomial order jump as well as the uniform solution are matching. The adaption does not incur any additional error. This is due to the element size being relatively large and in conjunction with the buffer of refined elements the coarse elements are situated far away from the interface region and thus the results match those of the uniform solution. What distinguishes the two solutions that make use of the adaption, is that the same accuracy can be achieved with fewer degrees of freedom for the case of the 2N/3 jump condition. Figures 4(a), 4(b) illustrate the field of the polynomial order for the case of Mesh 1 for two different time instants along with the 2N/3 criterion. Given a sufficient adaption rate the polynomial order jump does not affect the solution and does not cause any error for this specific mesh and simulation settings.   To understand the underlying difference between the two adapted solutions as well as the error generated during the polynomial adaption process, we show a cross section along the middle line of the bubble.

Spinodal Decomposition
The spinodal decomposition is a classic problem of the Cahn-Hilliard equation [3]. The main idea of this test is to introduce noise of different frequencies to an initial state of φ = 0 and allow it to evolve to the state of minimum energy. This is accompanied by a separation of the phases and the creation of bulk regions for each phase. We test the qualities and effectiveness of the adaption scheme on a two-and a three-dimensional spinodal decomposition.

Two-dimensional spinodal decomposition
The initial condition, φ 0 (x, y) =0.05 cos (0.105x) cos (0.11y) + (cos (0.13x) cos (0.087y)) 2 and the problem specification stem from the benchmark cases described in [75]. The visualization of the initial state of the simulation is presented in Figure 9(a). The parameters are specified in Table 6. The adaption frequency has been set to occur every 500 iterations. This problem has also been studied with the present numerical schemes in [8], and the sensitivity to the parameters of the time marching IMEX scheme has been identified. Also, the effect of having a structured or unstructured mesh has been quantified, and therefore the reader is invited to see the results therein. Here, we focus on the effect of the adaption using a more general unstructured mesh. The quantity of interest in the spinodal decomposition is the evolution of the free-energy. We also monitor the variation of the degrees of freedom as the simulation evolves.    The results presented in Figure 10(a) indicate that for that specific mesh the quantity of interest, the free-energy, converges when the polynomial order used is N ≥ 4. The final state is presented in Figure 9(b). Thus a polynomial order of N = 4 has been used as the basis for the comparison among the uniform and the adapted solution. The coarse level polynomial order used in this test case is N Coarse = 2. As presented in Figure 10(a), the adapted solution manages to achieve the same solution and the same final state. The case with the N − 1 criterion manages to achieve a 30% reduction in the degrees of freedom whereas the 2N/3 achieves a 35% reduction. A significant proportion of the simulation is carried out using fewer degrees of freedom.
The diagrams in Figure 11 present the field of the polynomial order as well as the position of the interfaces for different time instants of this spinodal decomposition problem. In Figure 11(a), which corresponds to t = 800, there are few elements which have been coarsened because the phases still undergo through the process of separation. As time evolves, the two phases are further separated as presented in Figure 11(b) and there is a significant decrease in the degrees of freedom. The final state is presented in Figure 11(c), which reveals that a large proportion of the elements has a coarse level polynomial order.
(a) t = 800 (b) t = 1600 (c) t = 10 4 Figure 11: visualization of the polynomial order distribution on the solution for (a) t = 800, (b) t = 1600 and (c) t = 10 4 for the spinodal decomposition test case with the 2N/3 criterion. The continuous black lines denote the locations of the interface where |φ| = 0.9.
In this case, we solve the Cahn-Hilliard equation using the chemical-free energy (9) and we also use an interface width parameter k which is four times lower than that of [8]. The latter has been chosen in order to simulate a system with a smaller interface and showcase the scheme's capabilities. The parameters of the simulation are presented in Table 7. We use a fine polynomial order of N Fine = 4, a coarse polynomial order of N Coarse = 2 and the 2N/3 jump criterion. The adaption frequency has been defined to occur at an interval of 1000 iterations.  The results presented in Figure 12 show that the results for the free-energy from the adapted solution match exactly those recovered when using a uniform polynomial order across the domain. As the solution evolves and the phases are separated, the adaption scheme reduces the degrees of freedom as presented in Figure 13. For this problem specification, the reduction achieved is 48%.
In Figure 14 the polynomial order distribution for different time instants is presented. The final state, as presented in Figure 14(d), is a flat interface separating the two phases.

Conclusions
In this work, we have established the theoretical framework for a free-energy stable discontinuous Galerkin scheme for the Cahn-Hilliard equation that targets p-non-conforming meshes. The analysis carried out has proven that for the continuous and the fully-discrete time setting, the stability characteristics of the conforming scheme of [8] can be maintained when transitioning to elements with non-uniform polynomial order. To do so, we use the BR1 numerical flux [68] and the made use of the standard mortar method [46]. In addition, we have designed a methodology for performing local polynomial refinement for the Cahn-Hilliard equation. This specific algorithm has the advantage of being simple and lowcost, based on the advantages that the higher order Discontinuous Galerkin Spectral Element Methods (DGSEM) offer as well as the Cahn-Hilliard model. We have also implemented and tested two different criteria to ensure that the scheme maintains the robustness and accuracy of the uniform solver. Lastly, the effectiveness of the method has been tested and verified through various two-and three-dimensional test cases, such as the spinodal decomposition, which have showcased that there is a significant reduction of the degrees of freedom required to attain the same accuracy as in the conforming version of the scheme.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.