A two-dimensional stabilized discontinuous Galerkin method on curvilinear embedded boundary grids
AA TWO-DIMENSIONAL STABILIZED DISCONTINUOUS GALERKIN METHOD ONCURVILINEAR EMBEDDED BOUNDARY GRIDS
ANDREW GIULIANI ∗ Abstract.
In this work, we propose a state redistribution method for high order discontinuous Galerkin methodson curvilinear embedded boundary grids. State redistribution relaxes the overly restrictive CFL condition that resultsfrom arbitrarily small cut cells when explicit time steppers are used. Thus, the scheme can take take time steps that areproportional to the size of cells in the background grid. The discontinuous Galerkin scheme is stabilized by postprocessingthe numerical solution after each stage or step of an explicit time stepping method. This is done by temporarily coarsening,or merging, the small cells into larger, possibly overlapping neighborhoods using a special weighted inner product. Then,the numerical solution on the neighborhoods is refined back onto the base grid in a conservative fashion. The advantageof this approach is that it uses only basic mesh information that is already available in many cut cell codes and does notrequire complex geometric manipulations. We prove that state redistribution is conservative and p -exact. Finally, we solvea number of test problems that demonstrate the encouraging potential of this technique for applications on curvilinearembedded geometries. Numerical experiments reveal that our scheme converges with order p + 1 in L and between p and p + 1 in L ∞ . Key words. state redistribution, discontinuous Galerkin methods, cut cell grids, explicit time stepping, curvilinearembedded boundaries.
AMS subject classifications.
1. Introduction.
Practical problems in computational fluid dynamics require the solution of hy-perbolic conservation laws on complex domains. Typically, the domain is first discretized into a mesh ofelements, then algorithms such as the finite volume (FV) or discontinuous Galerkin (DG) methods areused to approximate the exact solution on that mesh. When the boundary of the domain is complex,e.g., an aircraft, the most time-consuming aspect of the workflow can be the mesh generation phase [1].One technique that solves this problem is the use of embedded-boundary grids whereby the complexboundary is superimposed on a Cartesian grid. The computational mesh is then composed of Cartesiancells on the domain interior and irregular “cut” cells on the domain boundary. These mesh generationalgorithms are robust and automated.A problem that arises when using explicit time stepping methods on embedded boundary grids isthe small cell problem. Explicit methods require a time step that is proportional to the cell size inthe mesh. Since irregular cut cells on the domain boundary can be arbitrarily small, this results in anoverly stiff system of ordinary differential equations that would require impractically small time steps.A special treatment of these cut cells is therefore necessary for explicit time integrators to be used.Many solutions to the small cell problem in finite volume methods have been proposed, such asflux redistribution [2], h -box methods [3, 4], implicit/explicit time stepping [5, 6], dimensionally-splitflux stabilization [7, 8], cell merging [9, 10], and more recently state redistribution (SRD) [11]. Fluxredistribution is an approach where each cell is updated with a locally stable, but possibly nonconser-vative, time step [2]. Excess flux is then locally redistributed in order to maintain conservation. Thisis straightforward to implement in three dimensions, but is only first order accurate at the embeddedboundary. H -box methods are second order accurate, stable, and have desirable theoretical properties[3, 4]. This is achieved by increasing the domain of dependence of fluxes on cut cells, however, doing soin three dimensions seems difficult to implement. Implicit/explicit time integration schemes have alsobeen developed, where the small cells are integrated implicitly and the large cells are integrated explicitly[6]. Finally, cell merging is popular in two dimensions [9, 10]. However, similar to h -box methods, thistechnique suffers from implementation difficulties in three dimensions due to the many different, possiblyincompatible ways to create merging neighborhoods.Similar approaches have been applied to cut cell DG discretizations. Cell merging [12, 13, 14, 15],subcycling in time [13], and penalization terms [16] have been used to deal with the small cells. When ∗ Courant Institute, New York University, 251 Mercer St., New York, NY 10012 ([email protected])1 a r X i v : . [ m a t h . NA ] F e b trong anisotropic refinement is required, triangular cut cell approaches for the Navier-Stokes equationshave also been examined in [17, 18]. Our approach to the small cell problem is a modification of thestate redistribution method, which has been successfully applied to second order finite volume methodsin two dimensions [11]. The idea behind state redistribution is to temporarily merge cut cells intolarger, possibly overlapping neighborhoods in a postprocessing step applied after each explicit time step.For finite volume schemes, the spatial accuracy of state redistribution is increased by reconstructing apolynomial on each merging neighborhood. Then, these neighborhood polynomials are recombined in aconservative manner back onto the base grid. Applying state redistribution to DG numerical solutionsfollows a similar strategy, however, since the DG scheme does not rely on reconstructions, there arenotable differences. Specifically, the DG numerical solution is composed of discontinuous, locally definedpolynomials that are written as linear combinations of polynomial basis functions. Thus, we adaptour state redistribution approach to DG methods by associating to each possibly overlapping, mergingneighborhood a tailored polynomial basis. The DG solution is projected onto these overlapping mergingneighborhoods using a particular weighted inner product. Then, the numerical solution is projectedback onto basis functions defined on elements of the embedded boundary grid. The state redistributionmethod for finite volume methods is an attractive technique as it only requires mesh information thatis already available in many cut cell cell codes: the geometry of the cut cell, and its connectivityinformation. There are two additional requirements for DG methods: (1) the ability to compute innerproducts on the cut cells, and (2) basis functions on the cut cells and neighborhoods. In this work, wewill describe both the DG scheme on the base cut cell grid and how to relax the time step restrictionon cut cell grids using state redistribution. Typically, numerical schemes on embedded boundary gridssuffer from a loss of accuracy at the boundary, due to the irregularity of the grid and we will examinethis on a number of numerical examples.The structure of this article is as follows. In section 2 we briefly describe the spatial discretiza-tion used to approximate solutions to hyperbolic conservation laws. In section 3, we recall the stateredistribution method for finite volume schemes and outline how to extend it to arbitrarily high orderdiscontinuous Galerkin methods. This will involve the projections that we discussed previously. For easeof understanding and to reduce the necessary notation, this section uses a model problem in one spacedimension. In sections 4 and 5, we go into more detail for two-dimensional computations. In section6, we show that the approach is p -exact and also conservative. In section 7, we present a number ofnumerical experiments. In particular, we examine the error at cut cells, confirming the results in [12, 16]that the cut cell methods suffer some loss of accuracy at the boundary.
2. The discontinuous Galerkin method on cut cell meshes.
We aim to approximate solutionsto hyperbolic conservation laws(2.1) u t + ∇ · F ( u ) = 0 on Ω × (0 , T ) , using the discontinuous Galerkin method coupled with explicit time steppers on cut cell meshes. In theabove, Ω ⊂ R is the spatial domain, T is the final time, u : Ω × (0 , T ) → R m , m ∈ N is a vector ofconserved quantities and F is a flux function. The domain Ω is discretized with a mesh of elements K i,j , by superim-posing the boundary of ∂ Ω on a Cartesian grid and deleting the solid portions of the grid. The geometryin two dimensions is defined as a function that takes a point ( x, y ) and returns whether it’s inside oroutside the domain. In the volume, the mesh comprises regular Cartesian cells of dimension ∆ x by ∆ y .On the boundary, there are irregular polygons with both straight and curved edges that are determinedby the intersections of the boundary, ∂ Ω, with the Cartesian grid. The x and y coordinates of the curvededge are polynomials of degree q that interpolate ∂ Ω at q + 1 points. For example, the cut cell on theright in Figure 2.1 has one curved edge, where q = 2. As it is shown in [19], discontinuous Galerkinmethods require at least piecewise quadratic approximations to curved boundaries. Thus, we alwayspair the second order DG method ( p = 1) with a piecewise quadratic approximation of the boundary( q = 2). Otherwise, we let p = q when p >
1. Each cut cell is assumed to have at most two, possibly i,j +1 K i,j K i +1 ,j ∂ Ω ∂ Ω Fig. 2.1 . On the left, we plot a zoom of a curvilinear cut cell grid at the boundary, where whole cells are white andcut cells are light blue. The boundary of the domain ∂ Ω is the bold black line. On the right, we show the quadrature rulefor p = 2 of degree 4 on K i,j . There are N vol i,j + 1 = 15 quadrature points on this cell. The intersection of the cell with theembedded boundary are indicated by black circles ( • ) and the volume quadrature points of the cut cell are red circles ( • ). curved, irregular edges associated to the embedded boundary ( ∂ Ω). For example, the cut cells in Figure2.1 (highlighted in blue), only have one irregular edge. Cells that have two irregular edges lie on sharpcorners of the embedded boundary (Figure 7.3). We also do not allow split or tunneled cells for ease ofcode development. This restriction is not fundamental and we will need to implement these features inthree dimensions for complicated geometries. We have made our mesh generation code written in Pythonavailable at https://github.com/andrewgiuliani/PyGrid-2D. This code can be used to reproduce all thehigh order embedded boundary grids used in this work, as well as generate other embedded boundarygrids in two dimensions.
We use the modal discontinuous Galerkin method as the base schemeon this cut cell grid. Multiplying (2.1) by a test function v ( x, y ) ∈ H ( K i,j ) and integrating on cell K i,j ,we obtain (cid:90) K i,j u t v dx dy + (cid:90) K i,j ∇ · F ( u ) v dx dy = 0 ∀ v ∈ H ( K i,j ) . Integrating the second term by parts and using the divergence theorem, we have (cid:90) K i,j u t v dx dy + (cid:90) ∂K i,j v F ( u ) · n dl − (cid:90) K i,j F ( u ) · ∇ v dx dy = 0 ∀ v ∈ H ( K i,j ) . The solution u and test function v are then restricted to S p ( K i,j ), the space of polynomials of degree lessthan or equal to p on K i,j . On each cell, Cartesian and cut, we use a non-tensor product, orthonormalbasis for S p ( K i,j ) that has N p + 1 basis functions, where N p = ( p + 1)( p + 2) / −
1. Alternatively, atensor product basis [12] can be used. Tensor and non-tensor product bases have respectively ( p + 1) and ( p + 1)( p + 2) / U i,j ( x, y, t ) = N p (cid:88) k =0 c i,j,k ( t ) ϕ i,j,k ( x, y ) , here c i,j,k ∈ R m is the k th solution coefficient associated to basis function ϕ i,j,k on cell K i,j , wherethe ϕ i,j,k are orthonormal with respect to the inner product(2.3) (cid:104) f, g (cid:105) K i,j = 1 | K i,j | (cid:90) K i,j f g dx dy, where | K i,j | is the volume of cell i, j . The modal discontinuous Galerkin method is then given by thesystem of ordinary differential equations(2.4) ddt c i,j,k = − | K i,j | (cid:34)(cid:90) ∂K i,j ϕ i,j,k F ∗ ( U − i,j , U + i,j ) · n dl − (cid:90) K i,j F ( U i,j ) · ∇ ϕ i,j,k dx dy (cid:35) , for k = 0 , . . . , N p , where U − i,j and U + i,j are the numerical solutions corresponding to K i,j , and a neigh-boring cell, respectively, located along the element boundary ∂K i,j . Finally, n is an outward facing unitnormal, and F ∗ is a numerical flux function.The semidiscretization (2.4) is integrated in time using a ( p + 1)th order accurate Runge Kutta timestepping algorithm. Since the cut cell grid can contain arbitrarily small cut cells, (2.4) can be stiff, andbe subject to an arbitrarily small time step restriction. State redistribution will allow the use of explicittime stepping algorithms on cut cell grids with a time step that is proportional to the Cartesian cell size(2.5) ∆ t (cid:18) | a | ∆ x + | b | ∆ y (cid:19) ≤ p + 1 , where | a | , | b | are the maximum wave speeds in the x and y directions [20, 21].The DG method above requires an orthogonal polynomial basis on each polygonal cut cell. Wenumerically compute orthogonal basis functions and their gradients at quadrature points using the QRfactorization of a Vandermonde-like matrix (section 4.2). One issue with cut cell DG approaches that arises in-dependently of the small cell problem is the evaluation of inner products on arbitrary curved polygons.Quadrature rules are required to approximate volume integrals on elements (cid:90) K i,j f ( x, y ) g ( x, y ) dx dy ≈ N vol i,j (cid:88) q =0 w i,j,q f ( x i,j,q , y i,j,q ) g ( x i,j,q , y i,j,q ) , where N vol i,j + 1 is the number of volume quadrature points on K i,j , and (cid:80) N vol i,j q =0 w i,j,q = | K i,j | . Someapproaches are based on the divergence theorem [22], hierarchical moment fitting [23, 24], “speckling” ofquadrature points [17, 18], and others [25, 26]. We use the software described in [27, 28] for quadraturerules on polygons with curved edges, but other quadrature generation codes can be used as well. Byconstruction, this algorithm always results in quadrature nodes interior to the polygon with positiveweights. Another well-known approach is to triangulate the polygon and apply standard quadraturerules on each subtriangle [12]. Typically, this algorithm results in too many quadrature points andweights and does not generalize to three dimensions since not all polyhedra can be triangulated withoutintroducing vertices. Finally, quadrature by triangulation can also be difficult to robustly use withcurvilinear embedded boundaries. Initially, our code used this technique by triangulating the cut cells insuch a way that only one triangle had a curved edge. Then, we applied quadrature rules on both curvedand straight-edge triangles by mapping them to a canonical triangle. We did not find a satisfactorytechnique to guarantee an invertible mapping from curved to canonical triangles, and so abandonedobtaining quadrature rules by triangulation.On all elements, both whole and cut, we use quadrature rules of degree 2 p . Surface integrals areapproximated using standard Gauss-Legendre rules of degree 2 p + 1, with p + 1 points. (cid:57) K (cid:57) K (cid:57) K K K K h h αh h αh h h cell sizes ˆ K (cid:57) ˆ K ˆ K (cid:57) ˆ K (cid:57) ˆ K ˆ K ˆ K K (cid:57) K (cid:57) K (cid:57) K K K K Fig. 3.1 . Nonuniform grid for model problem used to describe the state redistribution method for both finite volumeand discontinuous Galerkin methods in Section 3. The top figure indicates the cell sizes. The bottom figure shows themerging neighborhoods and the overlap counts on each cell in the base grid. ˆ K i refers to the merging neighborhoodassociated to cell K i in the base grid. Neighborhoods that contain only one cell in the base grid are drawn in blue, whilethose that contain multiple, temporarily merged cells are drawn in red.
3. State redistribution in one dimension.
In this section, we demonstrate state redistributionwith a simple example that solves the linear advection equation(3.1) u t + au x = 0 , a > h i = h and on the small cells,and h i = αh for 0 < α <
1. Five cells (indexed by (cid:57) (cid:57)
2, 0, 2, 3) are large with size h and the remainingtwo cells (indexed by (cid:57) αh . For completeness, we first recall how to stabilizea first order upwind finite volume scheme using the approach described in [11]. Note that this schemecan be viewed as a discontinuous Galerkin method when p = 0. Then, we will describe how to extendthis approach to arbitrarily high order discontinuous Galerkin methods. State redistribution with finite volume methods.
We discretize (3.1) using the first order,upwind finite volume scheme(3.2) ˆ U i = U ni − a ∆ th i ( U ni − U ni − ) . To apply state redistribution to the above finite volume scheme, we must execute a simple geometricpreprocessing stage on the grid before the iterative portion of the finite volume solver.During this preprocessing stage, we associate a merging neighborhood ˆ K i to each cell in the basegrid, K i . Cells that are smaller than half the regular grid size ( h i < h/ K (cid:57) = K (cid:57) ∪ K (cid:57) ∪ K and ˆ K = K ∪ K ∪ K (Figure 3.1,in red). Cells that are larger than half the regular grid size ( h i ≥ h/
2) do not require merging and wehave ˆ K i = K i (Figure 3.1, in blue). Additionally, each cell in the base grid must count the number ofneighborhoods that overlap it, N i . For example, since K is a member of three neighborhoods: ˆ K (cid:57) , ˆ K ,and ˆ K , its overlap count is N = 3. We have provided the overlap counts for the cells on the base gridin Figure 3.1. fter the preprocessing stage has been completed, we apply state redistribution as a postprocessingstep after each forward Euler update of (3.2). On each neighborhood, we compute a special weightedsolution average. For example, on ˆ K (cid:57) the weighted solution average is(3.3) ˆ Q (cid:57) = 1 h/ αh + h/ (cid:18) h U (cid:57) + αh ˆ U (cid:57) + h U (cid:19) , where both cell size and solution average in the above sums are divided by their associated overlapcounts. Similarly, on ˆ K the weighted solution average isˆ Q = 1 h/ αh + h/ (cid:18) h U + αh ˆ U + h U (cid:19) . On neighborhoods that only contain one cell, the weighted solution average isˆ Q i = ˆ U i for i = (cid:57) , (cid:57) , , , . Finally, these weighted solution averages are recombined back onto the base grid. The final update ofa cell in the base grid is given by the average of the neighborhood values that overlap it. For example,since K is overlapped by ˆ K (cid:57) , ˆ K , and ˆ K , we have(3.4) U n +10 = 13 ( ˆ Q (cid:57) + ˆ Q + ˆ Q ) .K (cid:57) and K are overlapped by two neighborhoods, so their final updates are(3.5) U n +1 (cid:57) = 12 ( ˆ Q (cid:57) + ˆ Q (cid:57) ) and U n +12 = 12 ( ˆ Q + ˆ Q ) . On the cells overlapped by only one neighborhood, the final update is simply(3.6) U n +1 i = ˆ Q i for i = (cid:57) , (cid:57) , , . We prove in [11] that the above scheme in (3.4)-(3.6), is conservative. Additionally, we show that it canbe extended to second order FV methods by reconstructing a slope on each merging neighborhood andusing that slope to recombine the neighborhood averages back onto the base grid in a linearity preservingfashion. In the following section, we describe how to extend this idea to discontinuous Galerkin methods.
State redistribution with discontinuous Galerkin methods.
We discretize (3.1) in spaceusing an upwind RK-DG scheme where the degree of the polynomial approximation is p . Assuming anSSP-RK scheme, the first forward Euler update can be written(3.7) ˆ c i,k = c ni,k − a ∆ th i (cid:20) U ni ( x i +1 / ) ϕ i,k ( x i +1 / ) − U ni − ( x i − / ) ϕ i,k ( x i − / ) − (cid:90) K i U ni ϕ (cid:48) i,k dx (cid:21) , for k = 0 . . . p, where ˆ c i,k are the provisionally updated solution coefficients, U ni is the DG solution on cell K i at time t n and ˆ U i is its provisionally updated (possibly unstable) counterpart U ni ( x ) = N p (cid:88) k =0 c ni,k ϕ i,k ( x ) , ˆ U i ( x ) = N p (cid:88) k =0 ˆ c i,k ϕ i,k ( x ) . Additionally, the time step is ∆ t = ah/ (2 p + 1) on all cells, and ϕ i,k are the orthonormal basis functionsof S p ( K i ) with respect to the inner product in (2.3). − (a) (b) ˆ ϕ (cid:57) , ˆ ϕ (cid:57) , ˆ ϕ (cid:57) , ˆ ϕ (cid:57) , ˆ ϕ (cid:57) , K (cid:57) K (cid:57) K ˆ K (cid:57) K (cid:57) ˆ ϕ (cid:57) , ˆ ϕ (cid:57) , ˆ ϕ (cid:57) , ˆ ϕ (cid:57) , ˆ ϕ (cid:57) , Fig. 3.2 . In a), we plot the basis functions (solid lines) that are orthogonal with respect to the weighted inner product (cid:104)· , ·(cid:105) ˆ K (cid:57) in (3.8) . Although the basis functions in a) appear to be translations and rescaling of the standard Legendrepolynomials, a closer inspection reveals that they are not. We have plotted the normalized Legendre polynomials (dottedlines) on the same axes to illustrate the difference. In b), we plot the orthogonal basis functions with respect to (cid:104)· , ·(cid:105) K (cid:57) ,which are simply translations and rescaling of the standard Legendre polynomials. Note that functions in b) are plottedon element K (cid:57) , which is the center element in a). Similar to the SRD stabilized FV scheme described earlier, we must complete a preprocessing stage,whereby merging neighborhoods are associated to each cell in the base grid. For the illustration here,we use the same merging neighborhoods as in the finite volume example. The postprocessing stagedescribed below is more involved for the DG scheme due to the spatial accuracy required and comprisestwo projections:1. The first step projects the provisionally updated numerical solution from the base grid, ˆ U i , ontobasis functions defined on merging neighborhoods. They are called, ˆ ϕ i,k and they span S p ( ˆ K i ).2. The second step projects the numerical solution from the merging neighborhoods back onto thebasis functions defined on the base grid. They are called ϕ i,k , and they span S p ( K i ).We emphasize that the neighborhood basis functions are not generally the same basis functions as theones on the base grid. In this example, ϕ i,k and ˆ ϕ i,k are not identical, and are defined on differentdomains. Specifically, the basis functions on the base grid ϕ i,k are orthogonal with respect to theinner product in (2.3). In contrast, the basis functions defined on the merging neighborhoods ˆ ϕ i,k are orthogonal with respect to a particular weighted L inner product. For example, ˆ ϕ (cid:57) ,k is the k thorthonormal basis function of S p ( ˆ K (cid:57) ) with respect to the inner product(3.8) (cid:104) f, g (cid:105) ˆ K (cid:57) = 1 h/ αh + h/ (cid:18) (cid:90) K (cid:57) f g dx + (cid:90) K (cid:57) f g dx + 13 (cid:90) K f g dx (cid:19) . Note that in the above inner product, the cell sizes are divided by their associated overlap counts.Additionally, each integral is divided by the overlap count associated to its domain of integration: thefirst integral is divided by two since N (cid:57) = 2, the second integral is divided by one since N (cid:57) = 1,and the final integral is divided by three since N = 3. In Figure 3.2, we plot both ϕ (cid:57) ,k and ˆ ϕ (cid:57) ,k for k = 0 , , . . . , ϕ (cid:57) ,k and ˆ ϕ (cid:57) ,k are defined on different domains.Second, the basis functions ϕ (cid:57) ,k on the base grid are simply translations and rescalings of the Legendrepolynomials, as expected, while the neighborhood basis functions ˆ ϕ (cid:57) ,k are not. e are now ready to compute the projection ˆ Q (cid:57) ( x ) of the provisionally updated DG solutions, ˆ U (cid:57) ,ˆ U (cid:57) , ˆ U , onto the merging basis for ˆ K (cid:57) :(3.9) ˆ Q (cid:57) ( x ) = N p (cid:88) k =0 ˆ q i,k ˆ ϕ (cid:57) ,k . Due to the orthogonality of the basis functions, the coefficients ˆ q i,k are written(3.10) ˆ q (cid:57) ,k = 1 h/ αh + h/ (cid:18) (cid:90) K (cid:57) ˆ U (cid:57) ˆ ϕ (cid:57) ,k dx + (cid:90) K (cid:57) ˆ U (cid:57) ˆ ϕ (cid:57) ,k dx + 13 (cid:90) K ˆ U ˆ ϕ (cid:57) ,k dx (cid:19) for k = 0 . . . p, This operation can be viewed as a local coarsening of the provisionally updated DG solution from thecells K (cid:57) , K (cid:57) , K onto the merging neighborhood ˆ K (cid:57) . Due to the projection operation, ˆ Q (cid:57) ( x ) definedin (3.9) is the closest polynomial of degree p to the piecewise defined provisionally updated solution,ˆ U (cid:57) , ˆ U (cid:57) , and ˆ U on ˆ K (cid:57) , measured in the weighted norm (cid:113) (cid:104)· , ·(cid:105) ˆ K (cid:57) that follows from (3.8). Also, notethat formula (3.10) is a generalization of the finite volume formula in (3.3), and simplifies to (3.3) when p = 0 since the provisional solution updates ˆ U (cid:57) , ˆ U (cid:57) , ˆ U are constants for finite volume methods. Incontrast, they are polynomials for the DG method in (3.10).Each merging neighborhood has its own weighted inner product, which depends on how cells aremerged and their associated overlap counts. For example, ˆ ϕ ,k is the k th orthonormal basis function of S p ( ˆ K ) with respect to the weighted L inner product (cid:104) f, g (cid:105) ˆ K = 1 h/ αh + h/ (cid:18) (cid:90) K f g dx + (cid:90) K f g dx + 13 (cid:90) K f g dx (cid:19) , where again the cell sizes and integrals are divided by their associated overlap counts. Additionally, thesolution coefficients of the projected solution on ˆ K are writtenˆ q ,k = 1 h/ αh + h/ (cid:18) (cid:90) K ˆ U ˆ ϕ ,k dx + (cid:90) K ˆ U ˆ ϕ ,k dx + 13 (cid:90) K ˆ U ˆ ϕ ,k dx (cid:19) . The merging neighborhood associated to large cells contain only one cell of the base grid. The innerproduct one such cell, e.g., ˆ K is (cid:104) f, g (cid:105) ˆ K = 1 h / (cid:18) (cid:90) K f g dx (cid:19) = 1 h (cid:90) K f g dx, where again we divide the cell volume and integral by their associated overlap counts. However, noticethat the overlap counts cancel out and the inner product simplifies. We can make the more generalstatement that on the neighborhoods that contain only one cell of the base grid, we have (cid:104) f, g (cid:105) ˆ K i = (cid:104) f, g (cid:105) K i since (cid:104) f, g (cid:105) ˆ K i = 1 h i /N i (cid:18) N i (cid:90) K i f g dx (cid:19) = 1 h i (cid:90) K i f g dx = (cid:104) f, g (cid:105) K i , for i = (cid:57) , (cid:57) , , , . ue to cancellation of the overlap counts. Thus, the solution coefficients of the numerical solution onthe merging neighborhoods associated to large cells in the base grid areˆ q i,k = 1 h i (cid:18)(cid:90) ˆ K i ˆ U i ˆ ϕ i,k dx (cid:19) = ˆ c i,k since ˆ K i = K i and the neighborhood inner product is the same as the base grid inner product. Thisimplies that ˆ ϕ i,k = ϕ i,k for k = 0 . . . p and i = (cid:57) , (cid:57) , , ,
3. In other words, preprocessing is onlyrequired in the neighborhood of cut cells.The final step is to project the stabilized neighborhood polynomials from the merging neighborhoods,ˆ Q i , back onto the base grid. For a given cell, K i , in the base grid, we can average the stabilizedneighborhood polynomials that overlap it and then project that average onto the base grid basis functions ϕ i,k associated to K i . For example, K is overlapped by three merging neighborhoods: ˆ K (cid:57) , ˆ K , ˆ K .Therefore, the stabilized solution coefficients on K after one stage of the SSP-RK method are written c n +10 ,k = 1 h (cid:90) K
13 ( ˆ Q (cid:57) + ˆ Q + ˆ Q ) ϕ ,k dx, for k = 0 . . . p. The solution coefficients on cells overlapped by two merging neighborhoods, K (cid:57) and K , can be written(3.11) c n +1 (cid:57) ,k = 1 h (cid:57) (cid:90) K (cid:57)
12 ( ˆ Q (cid:57) + ˆ Q (cid:57) ) ϕ (cid:57) ,k dx, and c n +12 ,k = 1 h (cid:90) K
12 ( ˆ Q + ˆ Q ) ϕ ,k dx, for k = 0 . . . p . Finally, the solution coefficients on cells overlapped by one merging neighborhood is c n +1 i,k = 1 h i (cid:90) K i ˆ Q i ϕ i,k dx, for k = 0 . . . p, i = (cid:57) , (cid:57) , , . We will prove in section 6 that state redistribution is p -exact and conservative.We have just described how to apply state redistribution to a single forward Euler stage of an SSPRunge Kutta method. For non-SSP explicit Runge Kutta methods, state redistribution can be appliedto the intermediate and final solutions. In the next section, we will describe the details of the pre- andpostprocessing steps of the state redistribution algorithm for the discontinuous Galerkin method in twodimensions. One-dimensional convergence test: to conclude this section, we solve the linear advection equation(3.12) u t + u x = 0 on [ (cid:57) , ,u ( x,
0) = sin( πx ) , with periodic boundary conditions on the nonuniform grid in Figure 3.3c), which is similar to the grid inFigure 3.1. This grid has 2 N + 3 elements where 2 N + 1 cells have size h and the remaining two elementscan be arbitrarily small with size αh , where 0 < α < K is centered at x = 0 and h = 2 / (2 N + 1 + 2 α ) with α = 10 − . We solve (3.12) using the DG scheme in (3.7) stabilizedby state redistribution, where the time step is(3.13) ∆ t = 0 . p + 1 h on all cells. The L and L ∞ errors at the final time T = 1 are provided in Figure 3.3, where we observe theexpected p +1 rate of convergence in both norms. Even though there are two cells that are vastly smallerthan the ones used to determine the time step in (3.13), state redistribution allows us to explicitly timestep in a stable manner. The Python code that generated this convergence study is available at https://github.com/andrewgiuliani/PyDGSRD-1D. It can be used with other one-dimensional conservationlaws and non-uniform grids. − − − − − − N L e rr o r P1 (2 . . . . . (a) − − − − − − N L ∞ e rr o r P1 (2 . . . . . (b)(c) x = − h · · · h h αh hx = 0 K αh h · · · h hx = 1 N cells N cells Fig. 3.3 . In a) and b), we provide the L and L ∞ errors respectively, at the final time T = 1 on the grid shown inc). State redistribution does not affect the rate of convergence at the small cells since we observe the expected p + 1 rateof convergence in both the L and L ∞ norms.
4. Preprocessing in two dimensions.
In two dimensions, there are three main preprocessingoperations. First, merging neighborhoods must be associated to each cell in the cut cells (section 4.1).This operation relies only on mesh information already available in many cut cell codes: cell connectivityand the geometric vertices of cut cells. Then, basis functions must be generated on the base grid (section4.2), and on the merging neighborhoods (section 4.3). This operation relies on quadrature rules for theevaluation of inner products on arbitrary polygons (section 2.3). Since we have only considered staticembedded boundaries in this work, these preprocessing steps are only done once. If moving boundariesare considered, the mesh preprocessing would need to be done every time the boundary is modified.
Merging neighborhoods on two-dimensional cut cells grids aregenerated following the same procedure as described in [11]. Each cell in the base grid, K i,j , is associatedto a merging neighborhood ˆ K i,j . This merging neighborhood is formed by grouping irregular boundarycells with other nearby cells in the direction closest to the boundary normal. This is done until theneighborhood volume satisfies(4.1) (cid:88) ( r,s ) ∈ M i,j | K r,s | ≥
12 ∆ x ∆ y, where M i,j is a set containing the indices of each cell belonging to the merging neighborhood ˆ K i,j . Forexample, in Figures 4.1a) and b) M i,j = { ( i, j ) , ( i, j + 1) } and M i,j +1 = { ( i, j + 1) } . In (4.1), the chosenthreshold of ∆ x ∆ y/ K i,j associated to K i,j in Figure 4.1a). Theneighborhood is highlighted in green and is generated when K i,j is merged in the direction closest theboundary normal (upward) with K i,j +1 . In contrast, K i,j +1 does not need to merge with any neighboringcells, since it is already large enough and satisfies (4.1). Thus, its merging neighborhood is composed ofonly itself, ˆ K i,j +1 = K i,j +1 , and is highlighted in Figure 4.1b). Finally, notice that neighborhoods ˆ K i,j and ˆ K i,j +1 overlap on cell K i,j +1 . i,j (a) K i,j +1 (b) K i (cid:48) ,j (cid:48) (c) K i (cid:48) ,j (cid:48) (d) Fig. 4.1 . The merging neighborhood of cells K i,j and K i,j +1 are highlighted in a) and b), respectively. Theseneighborhoods overlap on cell K i,j +1 . The normal merging neighborhood does not always satisfy the volume constraint in (4.1) . In c), we show one such neighborhood, associated to K i (cid:48) ,j (cid:48) . In d), K i (cid:48) ,j (cid:48) is merged instead with cells on the × tile centered on K i (cid:48) ,j (cid:48) . This alternative approach to merging results in a neighborhood that is large enough. The portionsof the grid that are highlighted in blue are solid and do not belong to the computational grid. For certain boundary configurations, the normal merging neighborhood will not satisfy (4.1). Thisis shown in Figure 4.1c), for the neighborhood associated to K i (cid:48) ,j (cid:48) . In this case, one could merge K i (cid:48) ,j (cid:48) with cells located on the 3 × K i (cid:48) ,j (cid:48) (Figure 4.1d), also called the central mergingneighborhood. On all cells in the base grid, whole and cut, DG requires apolynomial basis with which to represent the numerical solution.We prefer to use a basis that is orthogonal with respect to the inner product(4.2) (cid:104) f, g (cid:105) K i,j = 1 | K i,j | N vol i,j (cid:88) q =0 w i,j,q f ( x i,j,q , y i,j,q ) g ( x i,j,q , y i,j,q ) , a discrete approximation to (2.3). In the above, w i,j,q , x i,j,q , y i,j,q are the q th volume quadrature weightsand points on K i,j (section 2.3). The inner product (4.2) requires the value of the basis functions at thevolume and surface quadrature points in order to evaluate the numerical solution in (2.2) for the volumeand surface integrals in (2.4). The gradient of the basis functions is also needed at the volume quadraturepoints for the volume integral in (2.4). This can be done by taking the monomial basis, and applyingthe modified Gram-Schmidt (MGS) algorithm [30, 31]. A coordinate transformation can also be appliedto the initial non-orthogonal monomial basis to reduce the effects of finite precision computations [30]. x, y ) i,j, min ( x, y ) i +2 ,j, min ∆ x i,j ∆ y i,j ∆ y i +2 ,j ∆ x i +2 ,j (a) η ξ . . . . . . . . (b) Fig. 4.2 . The above figures illustrate how the bounding box of a cut cell is used to generate basis functions. In a), weindicate the bottom left coordinates ( ◦ ) and dimensions of the bounding box for cells K i,j and K i +2 ,j used for the basisfunction calculation on the grid in Figure 4.1a) and b). In b), we plot the points ( ξ, η ) i +2 ,j,q for q = 0 , , . . . (indicatedby • ) at which the basis functions P k ( ξ, η ) are evaluated to populate the initial Vandermonde-like matrix V i +2 ,j in (4.4) for element K i +2 ,j . In this work, we describe another approach, where we determine these basis function values using the QR factorization of a matrix that we describe below.The matrix that we factorize is based on a Vandermonde-like matrix, populated with the valuesof initial basis functions that are not necessarily orthogonal with respect (4.2). This initial basis { P k ( ξ, η ) } k =0 ,...,N p is chosen to span S p ([0 , ) and is orthonormal with respect to the inner prod-uct (cid:104)· , ·(cid:105) [0 , . It is generated symbolically using the Gram-Schmidt procedure, using the monomialsordered as { , ξ, η, ξ , ξη, η , . . . } . The first three of these functions are P ( ξ, η ) = 1 , P ( ξ, η ) = 2 √ ξ − / , P ( ξ, η ) = 2 √ η − / . Then, we evaluate P k ( ξ i,j,q , η i,j,q ) at each volume quadrature point ( x i,j,q , y i,j,q ) on cell K i,j using themapping(4.3) ( ξ i,j,q , η i,j,q ) = (cid:18) x i,j,q − x i,j, min ∆ x i,j , y i,j,q − y i,j, min ∆ y i,j (cid:19) , where ( ξ i,j,q , η i,j,q ) is the volume quadrature point ( x i,j,q , y i,j,q ) mapped to [0 , , ∆ x i,j is the width ofthe cell i, j ’s bounding box in x direction, and x i,j, min is the minimum x component of all the verticesof K i,j . ∆ y i,j and y i,j, min are similarly defined in the y direction. These quantities are illustrated inFigure 4.2 for the mesh in Figure 4.1a) and b). Computing the basis function values at the quadraturepoints uses the QR factorization of a matrix in an oblique inner product [32], and requires three steps:1. Assemble the initial values of P k ( ξ i,j,q , η i,j,q ) into a Vandermonde-like matrix as follows(4.4) V i,j = P ( ξ i,j, , η i,j, ) P ( ξ i,j, , η i,j, ) . . . P N p ( ξ i,j, , η i,j, ) P ( ξ i,j, , η i,j, ) P ( ξ i,j, , η i,j, ) . . . P N p ( ξ i,j, , η i,j, ). . . P ( ξ i,j,N vol i,j , η i,j,N vol i,j ) P ( ξ i,j,N vol i,j , η i,j,N vol i,j ) . . . P N p ( ξ i,j,N vol i,j , η i,j,N vol i,j ) .
2. Multiply the rows of V i,j by the square root of the quadrature weights and using MATLAB we ompute the reduced QR factorization of the product W / i,j V i,j = Q i,j R i,j , where W i,j = diag( w i,j, , w i,j, , . . . , w i,j,N vol i,j ).3. Compute the product W (cid:57) / i,j Q i,j , which contains the orthogonal basis function values at thequadrature points W (cid:57) / i,j Q i,j = ϕ i,j, ( x i,j, , y i,j, ) ϕ i,j, ( x i,j, , y i,j, ) . . . ϕ i,j,N p ( x i,j, , y i,j, ) ϕ i,j, ( x i,j, , y i,j, ) ϕ i,j, ( x i,j, , y i,j, ) . . . ϕ i,j,N p ( x i,j, , y i,j, ). . . ϕ i,j, ( x i,j,N vol i,j , y i,j,N vol i,j ) ϕ i,j, ( x i,j,N vol i,j , y i,j,N vol i,j ) . . . ϕ i,j,N p ( x i,j,N vol i,j , y i,j,N vol i,j ) . The mass matrix computed from these basis functions is the identity matrix:( W (cid:57) / i,j Q i,j ) T W i,j ( W (cid:57) / i,j Q i,j ) = Q Ti,j Q i,j = I. The information stored in R i,j can be used to determine the gradient of the orthogonal basis functions.For example, the x and y components of the gradient of the orthonormal basis functions at the volumequadrature points are given by 1∆ x i,j ∂V i,j ∂ξ R (cid:57) i,j , and 1∆ y i,j ∂V i,j ∂η R (cid:57) i,j , respectively. In the above formulas, the 1 / ∆ x i,j and 1 / ∆ y i,j multipliers are due to the mapping (4.3).Similarly, the values of the basis functions at the surface quadrature points are given by S i,j R (cid:57) i,j , where S i,j is a Vandermonde-like matrix containing the initial basis functions { P k ( ξ, η ) } k =0 ,...,N p evaluated atthe Gauss-Legendre quadrature points on each edge of K i,j mapped to [0 , using (4.3). We also need to generate an orthogonal polyno-mial basis on each merging neighborhood that is orthogonal with respect to the weighted inner product(4.5) (cid:104) f, g (cid:105) ˆ K i,j = 1 | ˆ K i,j | (cid:88) ( r,s ) ∈ M i,j N r,s (cid:90) K r,s f g dx dy, where N r,s is the number of overlapping neighborhoods on cell ( r, s ), and | ˆ K i,j | = (cid:88) ( r,s ) ∈ M i,j | K r,s | N r,s , is the weighted volume of the merging neighborhood. The quadrature points at which the mergingneighborhood basis is computed are the same quadrature points as the ones used on the cells in basegrid (section 2.3). These basis function values are computed by adapting (4.3) to the bounding box ofthe merging neighborhood. That is, we evaluate P k ( ξ r,s,q , η r,s,q ) at each of the volume quadrature points( x r,s,q , y r,s,q ) on the merging neighborhood ˆ K i,j , ∀ ( r, s ) ∈ M i,j , using the mapping(4.6) ( ξ r,s,q , η r,s,q ) = (cid:18) x r,s,q − ˆ x i,j, min ∆ˆ x i,j , y r,s,q − ˆ y i,j, min ∆ˆ y i,j (cid:19) where ( ξ r,s,q , η r,s,q ) is again the volume quadrature point ( x r,s,q , y r,s,q ) mapped to [0 , , ∆ˆ x i,j is thewidth of merging neighborhood i, j ’s bounding box in the x direction, and ˆ x i,j, min is the minimum x ˆ x, ˆ y ) i,j, min (ˆ x, ˆ y ) i +2 ,j, min ∆ˆ x i,j ∆ˆ y i,j ∆ˆ y i +2 ,j ∆ˆ x i +2 ,j (a) η ξ . . . . . . . . (b) Fig. 4.3 . The above figures illustrate how the bounding box of a merging neighborhood is used to generate theneighborhood basis functions. In a), we indicate the bottom left coordinates ( ◦ ) and dimensions of the bounding box formerging neighborhoods ˆ K i,j and ˆ K i +2 ,j used for the merging basis function calculation on the grid on the top row ofFigure 4.1. In b), we plot the points ( ξ, η ) i +2 ,j,q for q = 0 , , . . . (indicated by • ) at which the basis functions P k ( ξ, η ) are evaluated to populate the initial Vandermonde-like matrix for merging neighborhood ˆ K i +2 ,j . Note that the quadraturepoints on each subelement of the merging neighborhoods are the same as the ones used in section 4.2. component of all the vertices of ˆ K i,j . ∆ˆ y i,j and ˆ y i,j, min are similarly defined in the y direction. Thesequantities are illustrated in Figure 4.3 for two merging neighborhoods on the mesh on the top row ofFigure 4.1. We can quickly precompute the merging basis function values using a discrete inner productand a QR factorization, as described in section 4.2. The discrete weighted inner product is given by(4.7) (cid:104) f, g (cid:105) ˆ K i,j = 1 | ˆ K i,j | (cid:88) ( r,s ) ∈ M i,j | K r,s | N r,s N vol i,j (cid:88) q =0 w r,s,q f ( r r,s,q , s r,s,q ) g ( r r,s,q , s r,s,q ) , where ( w r,s,q , x r,s,q , y r,s,q ) are the quadrature weights and points on each cell in the merging neighbor-hood ( r, s ) ∈ M i,j , and again (cid:80) q w r,s,q = | K r,s | and w r,s,q > q = 0 . . . N vol i,j . This basis will be usedin the projection operations that stabilize the numerical solution, described in the next section.
5. State redistribution in two dimensions.
We are now in a position to describe the stateredistribution method on two-dimensional cut cell grids. It is applied after each forward Euler step in ahigh order SSP-RK time stepper. For non SSP-RK time steppers, the SRD operator can be applied tothe intermediate and final solutions. The algorithm follows the operations described in one dimension(section 3), where the numerical solution is provisionally updated, then stabilized by two projectionoperations. The first projection moves the provisionally updated numerical solution from the basegrid onto the merging neighborhoods (section 3, step 1). The second projection moves the mergingneighborhood solution back onto the base grid in a conservative manner (section 3, step 2).To begin, we compute a forward Euler step or an intermediate solution of the RK scheme using thesame ∆ t on all cells, both whole and cut. For example, a forward Euler step for a scalar conservationlaw is written(5.1) ˆ c i,j,k = c ni,j,k − ∆ t | K i,j | (cid:34)(cid:90) ∂K i,j ϕ i,j,k F ∗ ( U − i,j , U + i,j ) · n dl − (cid:90) K i,j F ( U i,j ) · ∇ ϕ i,j,k dx dy (cid:35) , where ˆ c i,j,k for k = 0 , . . . N p , are the provisionally updated, possibly unstable degrees of freedom, and∆ t is the stable time step given by the CFL condition in (2.5). ext, we project the provisionally updated numerical solution onto the weighted merging basisprecomputed in section 4.3, using the discrete weighted inner product (4.7). The solution on the mergingbasis is ˆ Q i,j ( x ) = N p (cid:88) k =0 ˆ q i,j,k ˆ ϕ i,k ( x, y )where the solution coefficients are(5.2) ˆ q i,j,k = 1 | ˆ K i,j | (cid:88) ( r,s ) ∈ M i,j N r,s (cid:90) K r,s ˆ ϕ i,j,k ˆ U r,s dx dy, due to the orthogonality of ˆ ϕ i,j,k , andˆ U r,s ( x ) = N p (cid:88) k =0 ˆ c i,j,k ϕ i,j,k ( x, y ) . Finally, the numerical solution on the neighborhoods is projected back onto the base grid by aver-aging the overlapping neighborhood solutions:(5.3) c n +1 i,j,k = 1 | K i,j | (cid:88) ( r,s ) ∈ W i,j N i,j (cid:90) K i,j ϕ i,j,k ˆ Q r,s dx dy, where c n +1 i,j,k is the stabilized solution coefficient on the base grid, and W i,j is the set of neighborhoodindices that overlap cell ( i, j ). In the above formula, we average the projection from different, overlappingneighborhoods that contain the same cell in the base grid. As an example, the set of neighborhood indicesthat overlap K i,j +1 in Figure 4.1a) and b) is W i,j +1 = { ( i, j ) , ( i, j + 1) } . This set does not in fact needto be precomputed and the action of (5.3) can be implemented using a nested for loop (Algorithm 5.1). Note : the state redistribution method is a linear operation, so it can be implemented with a simplematrix vector product.
Algorithm 5.1
Implementation of (5.3) without constructing W i,j for i, j, k do c n +1 i,j,k ← end forfor i, j dofor ( r, s ) ∈ M i,j dofor k = 0 . . . N p do c n +1 r,s,k ← c n +1 r,s,k + ( (cid:82) K r,s ϕ r,s,k ˆ Q i,j dx dy ) / ( | K r,s | N r,s ) end forend forend for6. Accuracy and conservation. In this section, we show that state redistribution does not modifypolynomial solutions to (2.4), or in other words, it is p -exact. We also show that state redistribution isconservative. Similar to finite volume schemes [11], this DG scheme stabilized by state redistribution isnot monotone and as a result is not total variation diminishing. P -exactness follows from the fact that each step in the state re-distribution algorithm preserves polynomials of degree p . Consider the case where the provisionallyupdated DG solution is a polynomial of degree p on the entire domain ˆ U i,j ( x, y ) = f ( x, y ), with f ∈ S p (Ω), ∀ K i,j . Projecting the provisionally updated numerical solution onto each weighted merging asis using the weighted inner product (cid:104)· , ·(cid:105) ˆ K i,j , we have that the solution on all merging neighborhoodsis the original function ˆ Q i,j ( x, y ) = f ( x, y ), since the merging basis { ˆ ϕ i,j,k } k =0 ...N p spans S p ( ˆ K i,j ), ∀ K i,j . Projecting the average of overlapping neighborhood solutions back onto the base grid usingthe inner product (cid:104)· , ·(cid:105) K i,j , the solution in the final update (5.3) remains U n +1 i,j ( x, y ) = f ( x, y ), sincespan k =0 ...N p { ϕ i,j,k } = S p ( K i,j ), ∀ K i,j . Conservation follows for the same reason as for the finite volumescheme. We adapt the proof presented in [11] to DG methods. The total mass of the DG numericalsolution after one forward Euler step of the RK scheme is(6.1) M n +1 = (cid:88) i,j (cid:90) K i,j U n +1 i,j ( x, y ) dx dy = (cid:88) i,j | K i,j | c n +1 i,j, , because ϕ i,j, = 1 and orthogonality of the basis functions. Using the final update (5.3) for the solutioncoefficient c n +1 i,j, associated to the constant basis function in (6.1), we obtain(6.2) M n +1 = (cid:88) i,j N i,j (cid:88) ( r,s ) ∈ W i,j (cid:90) K i,j ˆ Q r,s dx dy. Note that the total mass can also be written as a sum of contributions from each merging neighborhoodby rearranging the sum in (6.2), i.e.,(6.3) M n +1 = (cid:88) i,j ˆ M i,j , where the contribution of merging neighborhood ( i, j ) is(6.4) ˆ M i,j = (cid:88) ( r,s ) ∈ M i,j N r,s (cid:90) K r,s ˆ Q i,j ( x, y ) dx dy, = | ˆ K i,j |(cid:104) ˆ Q i,j , (cid:105) ˆ K i,j . Recognizing that the first equation in (6.4) is the weighted projection of ˆ Q i,j ( x, y ) onto the neighbor-hood’s constant basis function ˆ ϕ i,j, = 1 scaled by | ˆ K i,j | , we have(6.5) ˆ M i,j = | ˆ K i,j | ˆ q i,j, , Using (5.2) in (6.5), we obtain(6.6) ˆ M i,j = (cid:88) ( r,s ) ∈ M i,j N r,s (cid:90) K r,s ˆ U r,s dx dy. Substituting the neighborhood contribution (6.6) into the expression for total mass on the grid (6.3), wehave(6.7) M n +1 = (cid:88) i,j (cid:88) ( r,s ) ∈ M i,j N r,s (cid:90) K r,s ˆ U r,s dx dy. This simplifies to(6.8) M n +1 = (cid:88) i,j (cid:90) K i,j ˆ U i,j dx dy, and shows that the mass is the same before and after state redistribution. . Numerical experiments. In this section, we demonstrate the performance of state redistri-bution method on a number of test problems. We solve a time-dependent scalar advection problem aswell as steady state problems in gas dynamics. In all numerical experiments, we use a time step that isproportional to the size of cells in the background mesh, given by∆ t (cid:18) | a | ∆ x + | b | ∆ y (cid:19) = 0 . p + 1 , where | a | , | b | are the maximum wave speeds in the x and y directions. On all embedded boundary grids,∆ x = ∆ y and N is the number of cells in each direction. We compute the discrete L and L ∞ errors(7.1) L = (cid:88) i,j N vol i,j (cid:88) q =0 w i,j,q | U i,j ( x i,j,q , y i,j,q ) − u ( x i,j,q , y i,j,q ) | ,L ∞ = max i,j max q | U i,j ( x i,j,q , y i,j,q ) − u ( x i,j,q , y i,j,q ) | , where U i,j ( x, y ) and u ( x, y ) are respectively the numerical and exact solutions. In this example, we solve for the solid body rotation of a pulse aroundan annulus described in [3, 33]. The conservation law u t − . π [( y − . u ] x + 0 . π [( x − . u ] y = 0 , rotates the initial condition u ( x, y,
0) = w ( θ − π/ , where w ( θ ) = 12 (cid:20) erf (5( π/ − θ )) + erf (5( π/ θ )) (cid:21) . The initial condition is plotted in Figure 7.2c) on the domain enclosed by two concentric discs of radii R = 0 .
75 and R = 1 .
25. We use the upwind numerical flux and impose zero flow at the boundary, i.e., F ∗ · n = 0 at each boundary quadrature point. Figure 7.2a) shows the computational domain embeddedon a 50 ×
50 grid. The exact solution is the initial condition rotated about the point (1 . , . T = 5 corresponds to one solid body rotation. The annulus is embedded on the domain[0 , . , where the x and y domain length is slightly perturbed from 3 to prevent cell degeneracies.The minimum volume fraction, min i,j | K i,j | / (∆ x ∆ y ), in this convergence test is 6.84e-10.In Figure 7.1d) and e), we plot the DG- P Q ×
10 grid at the final time. Segments corresponding to different elements areplotted in different colors. We see that there are no glitches or instabilities in the numerical solution onthe cut cells.
Discussion of the accuracy.
In Figures 7.2b) and d), we provide the L and L ∞ errors in (7.1),respectively. The expected p + 1 rate of convergence in the L norm is observed, but the rate ofconvergence in the L ∞ norm is between p and p + 1. This is in contrast to the one-dimensional results insection 3, where the full rate of convergence in both norms was observed. This is the price of automaticmesh generation that we are willing to pay. The nonsmoothness of the error in the L ∞ norm has beenobserved before and is due to the irregularity of the grid at the cut cells [11, 16]. In [16], it is found thatthe rate of convergence in L ∞ depends on the geometry of the domain. Here, we find that this rate alsodepends on the polynomial degree of approximation and varies between p and p + 1. Comparison with finite volume schemes.
We compare the numerical solution obtained withour DG scheme to the one obtained with a second order rotated grid h -box method for finite volumesdescribed in [33]. π/ π/ π/ π . . . θ (d) Outer bdry. π/ π/ π/ π . . . θ (e) Inner bdry. N FV ( h -box [33]) DG- P Q P Q (a) Domain E d N FV ( h -box [33]) DG- P Q P Q (b) Outer boundary E b N FV ( h -box [33]) DG- P Q P Q (c) Inner boundary E b Fig. 7.1 . On the left, we compare solutions to the solid body rotation problem in section 7.1 computed with DG-P1Q1, DG-P1Q2, and FV h-box [33] methods. In a), b), and c) the L errors on the domain, and on the inner andouter boundaries ( E d , E b ) of the numerical solution computed by all three numerical methods are provided, respectively.The errors in the FV column are taken from Table 3.1 of [33], and the numbers in parentheses are the estimated orderof convergence. On the right in d) and e), we plot the DG- P Q numerical solution on the outer and inner embeddedboundary of a × grid, respectively. Segments corresponding to different elements are plotted in different colors. In [33], the L error at the final time for second order accurate finite volume methods is measuredon the domain interior and boundary using(7.2) E d = (cid:80) i,j | U i,j − u ( x i,j , y i,j ) | | K i,j | (cid:80) i,j | u ( x i,j , y i,j ) | | K i,j | and E b = (cid:80) ( i,j ) ∈ B | U i,j − u ( x i,j , y i,j ) | | b i,j | (cid:80) ( i,j ) ∈ B | u ( x i,j , y i,j ) | | b i,j | , respectively, where ( x i,j , y i,j ) are the cell centroids, B is the set of indices of all cells that lie on thedomain boundary, and | b i,j | is length of the cut cell’s boundary segment. Since second order finite volumemethods typically represent curved boundaries with linear segments, we compute the L errors in (7.2)for both DG-P1Q1 and DG-P1Q2 methods.We find that the error is smaller on the outer boundary than on the inner boundary in Figures 7.1b)and c). This agrees with the results in [33], where this phenomenon is attributed to the outer boundaryhaving more cells to resolve the sloped regions of the pulse. As expected, the DG-P1Q2 scheme is moreaccurate than DG-P1Q1 scheme due to the higher fidelity representation of the boundary. Additionally,the DG-P1Q1 scheme is approximately an order of magnitude more accurate than the FV h -box schemein Figure 7.1a). The boundary errors of the DG schemes are smaller than the FV ones, although thedifference between FV and DG schemes here is not as pronounced. Finally, we observe that the FVscheme is second order accurate at the boundary, while the DG-P1Q1 and DG-P1Q2 schemes are not. a) − − − − − N L e rr o r P1Q2 (2 . . . . . (b)(c) − − − − − N L ∞ e rr o r P1Q2 (1 . . . . . (d) Fig. 7.2 . In a) we plot the annulus domain superimposed on the × background grid. In c) we plot the isolinesof exact solution at the initial and final time. In b) and d) we plot L and L ∞ errors of the solid body rotation problem,respectively. In the legend, we indicate the degree of polynomial approximation on the cells and the boundary segmentsby P and Q , respectively. The rate of convergence quoted in parentheses is computed by a least squares fit, plotted witha straight line. We observe the expected p + 1 rate of convergence in the L norm, and a rate of convergence between p and p + 1 in the L ∞ norm. Additionally, the error in the L ∞ norm is noisy, due to the irregularity of the grid at the cutcells. In this example, we solve for smooth transonic flow through a curved channel.This test problem is useful to evaluate the accuracy of high order methods with curved boundaries[34, 35, 36, 37, 38]. The exact solution at a given point ( x, y ) is computed by solving (cid:18) x − J − s x (cid:19) + ( y − s y ) = 14 ρ k , or the speed of sound c using the method of bisection, where J = 1 c + 13 c + 15 c −
12 ln (cid:18) c − c (cid:19) ,q = (cid:115) − c ) γ − , k = (cid:115) /q − ρ ( x − J/ − s x ) ,ρ = c / ( γ − ,γ = 1 .
4, and ρ , q are the density and speed of the gas, respectively. Additionally, ( s x , s y ) = (1 . ,
0) areshifts to center the channel on the domain. The pressure, and x, y components of the flow velocity aregiven by P = c ργ , u = (cid:112) q − v , and v = q k , respectively. The boundaries of the channel are embedded on [0 , . , and are given by the curves x ( k, q ) = 12 ρ (cid:18) q − k (cid:19) + J s x ,y ( k, q ) = 1 kqρ (cid:114) − q k + s y . The reflecting walls are described by k = 0 . k = 1 .
2, while the top boundary is described by q = 0 .
5. The domain and exact solution are shown in Figure 7.4a) and c). The exact solution is used asthe initial condition, and as the ghost states at the top and bottom boundaries. Roe’s numerical flux isused on the interior element interfaces, and the flux on the reflecting walls is given by(7.3) F ∗ · n = P n x P n y , where n = ( n x , n y ) is the outward facing normal and P is the internal pressure at the boundary at eachquadrature point. The solver is run until the stopping criterion(7.4) max i,j || c n +1 i,j − c ni,j || ∞ < (cid:57) is satisfied. The minimum volume fraction, min i,j | K i,j | / (∆ x ∆ y ), in this convergence test is 1.25e-09.This test case is interesting for two reasons. First, the flow smoothly transitions from subsonic tosupersonic and the Mach number reaches an approximate maximum of 1 .
42 near the bottom right cornerof the domain. Second, the embedded boundary presents sharp corners. The cut cells on which thesecorners lie may not have a unique normal merging neighborhood. For example, in Figure 7.3a), thecut cell on the top left corner of the domain can either be merged downward or to the right. Insteadof choosing one of the two possible normal merging neighborhoods, we use the 3 × L and L ∞ errors in the entropy in Figure 7.4b) and d), where the entropy is givenby P/ρ γ . We observe the expected p + 1 rate of convergence in the L norm and a rate between p and p + 1 in the L ∞ norm. This test problem was solved in [34], where their spectral volume scheme onunstructured triangles may form shocks that do not allow convergence of the method. For all meshesconsidered here, our solver converged and the stopping criterion (7.4) was attained. Finally, in [39] asecond order finite volume method solved a similar Ringleb flow problem on cut cell grids, and observeda rate of convergence of 2.02 in L and 1.40 in L ∞ . a) K i,j K i,j (b) K i − ,j − K i − ,j − (c) (a) K i,j K i,j (b) K i − ,j − K i − ,j − (c) Fig. 7.3 . In a), we plot the overlap counts of the Ringleb domain on the × grid. On most cut cells, we use thenormal merging neighborhood. When the volume constraint (4.1) cannot be satisfied, or there is a corner on the cut cell,we use the × central merging neighborhood. This is illustrated with the neighborhood associated to K i,j in b), whichoverlaps with the neighborhood associated to K i − ,j − shown in c). In both b) and c), the cells highlighted in green arethe merging neighborhoods and the portions that are highlighted in blue are solid and do not belong to the computationalgrid. The highlighted red cell in a), K i,j − , is overlapped by three neighborhoods: ˆ K i,j , ˆ K i − ,j − , and ˆ K i,j − . In this example, we compute the scattering ofa pressure pulse surrounded by five circular obstacles. The computational domain is [0 , and theobstacles are centered at (10 + 5 cos( ϕ k ) ,
10 + 5 sin( ϕ k )) for ϕ k = 2 πk/ − π/
3. The initial condition isgiven by ρuvP = − /γ + P /γ + 10 − exp( − b (cid:2) ( x − + ( y − (cid:3) ) , where γ = 1 . b = log(2) / . , and the disturbance in the pressure is a Gaussian located at (10 , ×
53 cut cell grid, where each boundary segment is a polynomialof degree five ( q = 5). On this grid, the minimum volume fraction, min i,j | K i,j | / (∆ x ∆ y ), is 7.28e-06.Roe’s numerical flux is used on the interior element interfaces, the flux along the reflecting obstacles isgiven by (7.3), and the ghost state on the outer boundaries of the domain is set to quiescent gas, ρ = 1, p = 1 /γ , and u = v = 0.We provide the DG- P Q T = 6 in Figure 7.5a) and c),respectively. The solution presents the same five-fold symmetry as the obstacles. The solution along theobstacle indicated by × is plotted in Figures 7.5b) and d), where we observe sub-cell resolution of thecomplex scattering phenomenon. a) − − − − − N L e rr o r P1Q2 (1 . . . . . (b)(c) − − − − − N L ∞ e rr o r P1Q2 (1 . . . . . (d) Fig. 7.4 . In a) we plot the Ringleb domain superimposed on a × background grid. In c) we plot the Machisolines of the exact solution, where the sonic line is labeled. In b) and d) we plot L and L ∞ entropy errors, respectively.In the legend, we indicate the degree of polynomial approximation on the cells and the boundary segments by P and Q ,respectively. The rate of convergence quoted in parentheses is computed by a least squares fit, plotted with a straight line.
8. Conclusions.
We have presented a state redistribution method for discontinuous Galerkin meth-ods on curvilinear embedded boundary grids. This technique preserves the formal order of accuracy ofthe base DG scheme and is conservative. The advantage of state redistribution is that it only uses basicmesh information that is already available in many cut cell codes and does not require complex geometriccomputations. Numerical examples reveal that the order of accuracy of the DG scheme is p + 1 in the L norm and between p and p + 1 in L ∞ . We find that rate in L ∞ depends on the problem, the geometryof the embedded boundary, as well as the polynomial degree of approximation p . This is a price thatwe are willing to pay in exchange for automatic mesh generation. Future work includes studying moreclosely the loss of accuracy at the boundary and attempt to improve it. We will also extend the methodto three dimensions and more realistic engineering geometries. In particular, we will have to considerseverely ill-shapen cut cells that commonly occur in three dimensions. Finally, when we extend ourmethod to viscous flows, we would like to adapt state redistribution to the Navier-Stokes equations.
9. Code availability.
We have made a one-dimensional DG-SRD implementation in Python avail-able at https://github.com/andrewgiuliani/PyDGSRD-1D. This code will reproduce the convergence a) − π − π/ π/ π − − · − θ ( P − / γ ) / − (b) × Fig. 7.5 . In a), we plot the scaled deviation from the mean pressure, ( P − /γ ) / − , at the final time. In b), weprovide this pressure deviation along the obstacle indicated by a × . Segments corresponding to different elements areplotted in different colors. We observe that the numerical solution retains the five-fold symmetry of the domain. On theembedded boundary, we observe the complex variation of the pressure in the scattered pulse. test in section 3 and can be applied to other conservation laws on nonuniform grids. We have also madeour cut cell mesh generation code written in Python available at https://github.com/andrewgiuliani/PyGrid-2D. This code can be used to reproduce all the grids used in the numerical experiments insection 7, as well as generate other high order embedded boundary grids in two dimensions.
10. Acknowledgements.
The author would like to thank Marsha Berger, Sandra May, GeorgStadler, and Lilia Krivodonova for their reading of the manuscript and helpful comments. We wouldalso like to thank Marco Vianello for the helpful discussions on using the quadrature rule generationcode. AG is partially supported by the Natural Sciences and Engineering Research Council of Canada(NSERC) under Award No. PDF-546085-2020 and by the Simons Foundation/SFARI (560651, AB).
REFERENCES[1] J. Dongarra, J. Hittinger, J. Bell, L. Chacon, R. Falgout, M. Heroux, P. Hovland, E. Ng, C. Webster, and S. Wild,“Applied mathematics research for exascale computing,” tech. rep., Lawrence Livermore National Lab, 2 2014.[2] P. Colella, D. T. Graves, B. J. Keen, and D. Modiano, “A Cartesian grid embedded boundary method for hyperbolicconservation laws,”
Journal of Computational Physics , vol. 211, no. 1, pp. 347–366, 2006.[3] M. Berger and C. Helzel, “A simplified h-box method for embedded boundary grids,”
SIAM J. Sci. Comput. , vol. 34,2012.[4] C. Helzel, M. Berger, and R. LeVeque, “A high-resolution rotated grid method for conservation laws with embeddedgeometries,”
SIAM J. Sci. Comput. , vol. 26, pp. 785–809, 2005.[5] S. Jebens, O. Knoth, and R. Weiner, “Linearly implicit peer methods for the compressible Euler equations,”
AppliedNumerical Mathematics , vol. 62, no. 10, pp. 1380 – 1392, 2012. Selected Papers from NUMDIFF-12.[6] S. May and M. Berger, “An explicit implicit scheme for cut cells in embedded boundary meshes,”
J. Sci. Comput. ,2016.[7] R. Klein, K. Bates, and N. Nikiforakis, “Well-balanced compressible cut-cell simulation of atmospheric flow,”
Philo-sophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences , vol. 367,no. 1907, pp. 4559–4575, 2009.[8] N. Gokhale, N. Nikiforakis, and R. Klein, “A dimensionally split Cartesian cut cell method for hyperbolic conservationlaws,”
Journal of Computational Physics , vol. 364, pp. 186–208, 2018.[9] J. J. Quirk, “An alternative to unstructured grids for computing gas dynamic flows around arbitrarily complextwo-dimensional bodies,”
Computers & fluids , vol. 23, no. 1, pp. 125–142, 1994.[10] S. Xu, T. Aslam, and D. Stewart, “High resolution numerical simulation of ideal and non-ideal compressible reacting23ows with embedded internal boundaries,”
Combustion Theory and Modelling , vol. 1, no. 1, pp. 113–142, 1997.[11] M. Berger and A. Giuliani, “A state redistribution algorithm for finite volume schemes on cut cell meshes,”
Journalof Computational Physics , vol. 428, p. 109820, 2021.[12] R. Qin and L. Krivodonova, “A discontinuous Galerkin method for solutions of the Euler equations on Cartesian gridswith embedded geometries,”
Journal of Computational Science , vol. 4, no. 1, pp. 24 – 35, 2013. ComputationalMethods for Hyperbolic Problems.[13] S. Schoeder, S. Sticko, G. Kreiss, and M. Kronbichler, “High-order cut discontinuous Galerkin methods with localtime stepping for acoustics,”
International Journal for Numerical Methods in Engineering , vol. 121, no. 13,pp. 2979–3003, 2020.[14] D. Krause and F. Kummer, “An incompressible immersed boundary solver for moving body flows using a cut celldiscontinuous Galerkin method,”
Computers & Fluids , vol. 153, pp. 118–129, 2017.[15] B. M¨uller, S. Kr¨amer-Eis, F. Kummer, and M. Oberlack, “A high-order discontinuous Galerkin method for compress-ible flows with immersed boundaries,”
International Journal for Numerical Methods in Engineering , vol. 110,no. 1, pp. 3–30, 2017.[16] C. Engwer, S. May, A. N¨ußing, and F. Streitb¨urger, “A stabilized DG cut cell method for discretizing the lineartransport equation,”
SIAM Journal on Scientific Computing , vol. 42, no. 6, pp. A3677–A3703, 2020.[17] K. J. Fidkowski and D. L. Darmofal, “A triangular cut-cell adaptive method for high-order discretizations of thecompressible Navier–Stokes equations,”
Journal of Computational Physics , vol. 225, no. 2, pp. 1653 – 1672, 2007.[18] K. Fidkowski and D. Darmofal, “An adaptive simplex cut-cell method for discontinuous Galerkin discretizations ofthe Navier-Stokes equations,” in , p. 3941, 2007.[19] F. Bassi and S. Rebay, “High-order accurate discontinuous finite element solution of the 2D Euler equations,”
Journalof Computational Physics , vol. 138, no. 2, pp. 251–285, 1997.[20] B. Cockburn and C.-W. Shu, “The Runge–Kutta discontinuous Galerkin method for conservation laws V: Multidi-mensional systems,”
Journal of Computational Physics , vol. 141, no. 2, pp. 199 – 224, 1998.[21] R. Qin and L. Krivodonova, “Spectrum of the discontinuous Galerkin spatial discretization for two dimensionalproblems,”
Preprint , 2020.[22] Y. Sudhakar, J. M. De Almeida, and W. A. Wall, “An accurate, robust, and easy-to-implement method for integra-tion over arbitrary polyhedra: application to embedded interface methods,”
Journal of Computational Physics ,vol. 273, pp. 393–415, 2014.[23] B. M¨uller, F. Kummer, and M. Oberlack, “Highly accurate surface and volume integration on implicit domains bymeans of moment-fitting,”
International Journal for Numerical Methods in Engineering , vol. 96, no. 8, pp. 512–528, 2013.[24] F. Kummer, “Extended discontinuous Galerkin methods for two-phase flows: the spatial discretization,”
InternationalJournal for Numerical Methods in Engineering , vol. 109, no. 2, pp. 259–289, 2017.[25] S. Mousavi, H. Xiao, and N. Sukumar, “Generalized Gaussian quadrature rules on arbitrary polygons,”
InternationalJournal for Numerical Methods in Engineering , vol. 82, no. 1, pp. 99–113, 2010.[26] R. Saye, “High-order quadrature methods for implicitly defined surfaces and volumes in hyperrectangles,”
SIAMJournal on Scientific Computing , vol. 37, no. 2, pp. A993–A1019, 2015.[27] A. Sommariva and M. Vianello, “Computing Tchakaloff-like cubature rules on spline curvilinear polygons,”
DolomitesResearch Notes on Approximation ∼ alvise/software.html.[29] M. Berger, “A note on the stability of cut cells and cell merging,” Applied Numerical Mathematics , vol. 96, pp. 180– 186, 2015.[30] F. Bassi, L. Botti, A. Colombo, D. A. Di Pietro, and P. Tesini, “On the flexibility of agglomeration based physicalspace discontinuous Galerkin discretizations,”
Journal of Computational Physics , vol. 231, no. 1, pp. 45–65,2012.[31] L. Mascotto, “Ill-conditioning in the virtual element method: Stabilizations and bases,”
Numerical Methods forPartial Differential Equations , vol. 34, no. 4, pp. 1258–1281, 2018.[32] B. R. Lowery and J. Langou, “Stability analysis of QR factorization in an oblique inner product,” arXiv preprintarXiv:1401.5171 , 2014.[33] C. Helzel, M. J. Berger, and R. J. Leveque, “A high-resolution rotated grid method for conservation laws withembedded geometries,”
SIAM Journal on Scientific Computing , vol. 26, no. 3, pp. 785–809, 2005.[34] Z. Wang and Y. Liu, “Extension of the spectral volume method to high-order boundary representation,”
Journal ofComputational Physics , vol. 211, no. 1, pp. 154 – 178, 2006.[35] M. Dumbser, M. K¨aser, V. A. Titarev, and E. F. Toro, “Quadrature-free non-oscillatory finite volume schemeson unstructured meshes for nonlinear hyperbolic systems,”
Journal of Computational Physics , vol. 226, no. 1,pp. 204 – 243, 2007.[36] P. Houston and N. Sime, “Automatic symbolic computation for discontinuous Galerkin finite element methods,”
SIAM Journal on Scientific Computing , vol. 40, no. 3, pp. C327–C357, 2018.[37] Z. J. Wang, K. Fidkowski, R. Abgrall, F. Bassi, D. Caraeni, A. Cary, H. Deconinck, R. Hartmann, K. Hillewaert,H. T. Huynh, et al. , “High-order CFD methods: current status and perspective,”
International Journal forNumerical Methods in Fluids , vol. 72, no. 8, pp. 811–845, 2013.[38] L. Ivan and C. P. Groth, “High-order solution-adaptive central essentially non-oscillatory (CENO) method for viscous24ows,”
Journal of Computational Physics , vol. 257, pp. 830 – 862, 2014.[39] W. J. Coirier and K. G. Powell, “An accuracy assessment of Cartesian-mesh approaches for the Euler equations,”