A high-order conservative remap for discontinuous Galerkin schemes on curvilinear polygonal meshes

Journal of Computational Physics 399 (2019) 108931

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

A high-order conservative remap for discontinuous Galerkin schemes on curvilinear polygonal meshes Konstantin Lipnikov a,∗ , Nathaniel Morgan b a

Applied Mathematics and Plasma Physics Group, Theoretical Division, Mail Stop B284, Los Alamos National Laboratory, Los Alamos, NM 87545, United States of America b Methods and Algorithms Group, X-Computational Physics Division, Mail Stop B284, Los Alamos National Laboratory, Los Alamos, NM 87545, United States of America

a r t i c l e

i n f o

Article history: Received 19 December 2018 Received in revised form 14 June 2019 Accepted 4 September 2019 Available online 10 September 2019 Keywords: Advection-based remap High-order schemes Polygonal meshes Discontinuous Galerkin Curvilinear meshes Data remap

a b s t r a c t A data transfer (called later remap) of physical fields between two meshes is an important step of arbitrary Lagrangian-Eulerian (ALE) simulations. This step is challenging for highorder discontinuous Galerkin schemes since the Lagrangian flow motion leads to high-order meshes with curved faces. It becomes even more challenging for unstructured polygonal meshes that do not have a polynomial map from the reference to a current cell. We propose and analyze a new framework to create remap schemes on curvilinear polygonal meshes based on the theory of virtual element projectors. We derive a conservative remap scheme that is high-order accurate in space and time. The properties of this scheme are studied numerically for smooth and discontinuous fields on unstructured quadrilateral and polygonal meshes. Published by Elsevier Inc.

1. Introduction One way to increase the predictive power of Lagrangian hydrodynamic simulations, without an excessive mesh refinement, is to use high-order numerical schemes that work on unstructured polytopal (polygonal in 2D, polyhedral in 3D) meshes. In the arbitrary Lagrangian-Eulerian (ALE) framework, the Lagrangian step of a simulation is combined with a mesh optimization step followed by a data transfer step from the Lagrangian mesh to the modified (optimized) mesh [11,10,38]. The last two steps are as critical for the overall simulation accuracy as the first step, since they must preserve characteristic mesh features as well as important mathematical and physical properties of the Lagrangian solution, see e.g. [35,24,18]. It is appropriate to mention here that there exist other ALE frameworks that solve PDEs written in a moving frame other than the Lagrangian frame [21,23]. Different mesh velocity implies that these frameworks perform continuous mesh modification and have no explicit remap step. The focus of this paper is on algorithms for data transfer for numerical schemes using discontinuous Galerkin (DG) approximations. A straightforward approach would be to employ the L 2 finite element projectors. This approach is needed in the reconnection ALE [33] approach due to local change of the mesh topology. Their practical implementation requires the intersection of two polygonal meshes which is feasible in 2D provided that both meshes have straight edges. Extension of this algorithm to 3D can be done under the assumption of flatness of mesh faces, although its computational cost is much

*

Corresponding author. E-mail addresses: [email protected] (K. Lipnikov), [email protected] (N. Morgan).

https://doi.org/10.1016/j.jcp.2019.108931 0021-9991/Published by Elsevier Inc.

2

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

higher. The mesh-mesh intersection approach becomes impractical when cells in the Lagrangian mesh and/or the optimized mesh have curved faces which is one of the interesting features of high-order schemes for Lagrangian hydrodynamics. The alternative approach is to express the data transfer step as a dynamic process governed by a partial differential equation, more precisely, a linear transport equation. If mesh optimization makes small perturbations to the original mesh, on the order of the mesh size, then the remap may be written as a simple flux-form transport algorithm [17]. The swept region algorithms work under similar assumption [24]. We do not make such an assumption here. The PDE-type approach is natural for meshes with strongly curved faces; however, it requires one to define the advection velocity (we also call it the remap velocity) everywhere in the computational domain which is straightforward only in the case of finite element meshes [16,3]. Indeed, there exist explicit formulas for maps between reference and curved triangles or between a square and a quadrilateral. However, there exists no polynomial map between two arbitrary quadrilaterals. In the case of polygonal meshes, such maps can be developed by breaking mesh cells into triangular pieces and using quadrature formulas on triangles. Besides higher computational cost, this is not a DG approach in spirit. For hyperbolic PDEs, the first DG method was introduced by Reed and Hill in 1973 [39] to simulate neutron transport. The DG methods are typically locally conservative, stable, and high-order accurate methods that can handle complex geometries, polytopal meshes and polynomial approximations of different order in neighboring mesh cells. These properties have brought DG methods into a huge variety of applications including but not limited to compressible [7,15,31] and incompressible [37] fluid dynamics, granular flows [20], semiconductor device simulation [12], viscoelasticity [5], transport of a contaminant in porous media [1], and Lagrangian hydrodynamics on linear [34,43,29,31] and curvilinear [44,32] meshes. Some original applications of DG methods include the Hamilton-Jacobi [27] and Korteweg-deVries [40] equations. Application of the DG method to the level set problems is discussed e.g. in [22,30]. There exist a number of DG schemes for the transport equation, exploiting various numerical technologies, that could be used for the remap purpose. In many of these schemes the velocity is assumed to be given and the focus is on accuracy and monotonicity. A numerical scheme that preserves positivity of density and pressure in DG schemes for compressible Euler equations is described in [46]. The idea is to control the polynomial values at a set of quadrature points. It is not obvious how to extend this idea to polygonal meshes. Finally, a flux-corrected transport (FCT) maximum principle preserving scheme in [2] is applicable to both advection remap and transport problems. This scheme employs positivity preserving Bernstein polynomials, hence it is also limited to finite-element meshes. An interesting particle-based remap scheme is proposed in [34], where each particle carries mass, momentum and energy. This remap scheme is locally conservative and its accuracy can be adjusted, but its extension to polygonal meshes is not obvious at all. A DG scheme of order k + 1 uses polynomials of order k to represent physical fields. We assume that the exact remap velocity is a specially designed virtual function. The trace of this function is a vector polynomial of order k on each (possibly curved) mesh edge. Inside mesh cells, pointwise values of this function are not available. We propose to approximate the unknown remap velocity inside a mesh cell using a vector polynomial of order k + 1. One additional order allows us to recover exact maps between two cells in a few special cases. For instance, a bilinear map is needed to map a square to a quadrilateral. A biquadratic map is needed to map a square to a quadrilateral with edges described by parabolas. The high-order approximation of the remap velocity in a mesh cell is based on the theory of virtual element projectors. The virtual element method (VEM) is a non-trivial generalization of the finite element to polytopal meshes [8,4] which emerged recently from the mimetic finite difference method [9]. The projector produces a polynomial with the minimal energy norm which could be considered as a measure of its optimality. Another useful property of the VEM projector is that it remains numerically stable for non-convex and degenerate cells as long as these cells are shape regular. We present the new technology for the case of polygonal meshes with straight and curved edges. The numerical experiments indicate significant error reduction between the finite volume (the DG scheme of order one) and the higher order DG schemes even on meshes with moderate resolution. Data transfer must satisfy various physical and mathematical requirements. The list of requirements includes but is not limited to conservation laws, monotonicity preservation and local and global accuracy. Our weak formulation is written in a flux form, which guarantees the conservation law. Limiting solution towards its mean value is often needed to avoid the creation of artificial extrema. One interesting feature of the remap problem is that the Barth-Jespersen limiter [6] for highorder schemes violates the conservation law. This is a huge difference with the DG schemes for Lagrangian hydrodynamics, where this limiter does not affect the conservation law. We propose a post-processing algorithm that leads to a conservative scheme. The paper outline is as follows. In Section 2, we formulate the PDE model for remapping conservative fields between two meshes. In Section 3, we derive its weak formulation and present semi-discrete equations in reference coordinates. In Section 4, we give description of the fully discrete scheme and study its properties. In Section 5, we confirm our finding with a few numerical experiments on square and polygonal meshes. 2. Advection based remap We start with a brief derivation of a PDE model for the advection-based remap. Consider a computational (Lagrangian) ˜ h . We assume that both meshes cover the same polygonal domain  without mesh h and a modified (optimized) mesh  gaps and overlaps. In addition, we assume that these meshes have the same number of polygonal cells and each cell in h ˜ h with the same number of vertices. corresponds to a cell in 

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

3

Fig. 1. Illustration of a Lagrangian (left) and modified (right) meshes with five cells and some basic notations. Mesh vertices are marked with solid disks.

˜ h , respectively. Hereafter, we refer Let ξ = (ξ1 , ξ2 ) and x˜ = (˜x1 , x˜ 2 ) denote the position vectors on meshes h and  to ξ as the reference coordinate system. We introduce a bijective map x˜ = F(ξ ) on  that maps nodes and edges of h ˜ h , see Fig. 1. The mesh cells can be generalized polygons with curved edges including to respective nodes and edges of  non-convex and degenerate polygons. We make the following assumption. For a DG scheme of order k + 1, the map F(ξ ) on each mesh face f is a vector polynomial of order k.

(A1)

The bijective map induces a space-time map

x(τ , ξ ) = ξ + τ (F(ξ ) − ξ ), with dimensionless pseudo-time

u=

0 ≤ τ ≤ 1,

τ . The remap velocity u is defined naturally as the time derivative of this map:

∂x = F(ξ ) − ξ . ∂τ

(1)

Note that the remap velocity is independent of the pseudo-time. Similar to the physical motion, velocity u defines a particle trajectory. A change of field ρ along this trajectory is given by the material derivative. Since the field is not changed from the viewpoint of a fixed observer, we have

dρ dτ

=

∂ρ + u · ∇ρ = u · ∇ρ. ∂τ

(2)

Similar equations can be derived for vector fields. Thus, it is sufficient to analyze equation (2) for a scalar field. 3. Semi-discrete formulation Due to the one-to-one correspondence between mesh objects, we use generic symbol c (resp., f ) to denote a polygonal ˜ h . We refer to f as a face to stress the fact that the presented technology can be extended cell (resp., face) in either h or  ˜ h . Let n f (x) be to 3D. When the space-time position of a mesh cell is important, we write c (τ ), with c (0) ∈ h and c (1) ∈  the unit normal vector to face f at point x whose orientation is fixed once and for all. When the space-time orientation of the normal is clear from the context, we write n f . To shorten notations, we introduce N f = n f (ξ ). Finally, let |c (τ )| denotes the cell volume. In cell c, we define a basis ψic (x) which has simple polynomial representation at time τ = 0, see formula below. Let ξ c denote the centroid of cell c ∈ h . For the second-order DG scheme, the Taylor basis [41,42] includes three functions:

ψ0c (ξ ) = 1,

ψ1c (ξ ) = ac1 (ξ1 − ξc,1 ),

ψ2c (ξ ) = ac2 (ξ2 − ξc,2 ).

Higher-order schemes add to this basis high-order polynomials of the form

  ψic (ξ ) = aci (ξ1 − ξc,1 )αi (ξ2 − ξc,2 )βi − bci ψ0c ,

i > 2,

where αi and βi are non-negative integers. The scaling factors aci and the orthonormalization factors bci are calculated from these two conditions:





ψ0c c ( 0)

ψic

ψic ψic dV = |c (0)|,

dV = δ0,i , c ( 0)

where δ0,i is the Kronecker symbol. The scaling factors improve spectral properties of mass matrices that appear later in the discrete scheme. Hence, they increase the stability range of an explicit time integration scheme. We assume that the basis functions remain constant along particle trajectories, i.e.

4

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

dψic dτ

= 0.

(3)

A discrete space-time field

ρh in cell c has the following expansion in the Taylor basis:

ρh |c (x) = ρ0c + ρ1c ψ1c (x) + ρ2c ψ2c (x) + ρ3c ψ3c (x) + . . . . A DG scheme of order k + 1 uses polynomials of order k; hence, the above expansion has (k + 1)(k + 2)/2 terms. For a face f shared by two cells c 1 and c 2 and the normal vector pointing from c 1 to c 2 , we define the jump and average operators as follows:

[ψ] = ψ|c1 − ψ|c2 ,

{ψ} =

1 2

 ψ|c1 + ψ|c2 .

Let us multiply (2) by a test function ψ , integrate the right-hand side by parts, and apply the Reynolds transport theorem to the left-hand side. Since ρ and ψ could be discontinuous functions, we have to first perform calculations cell-by-cell and then sum up the results:



d





ρ ψ dv = −

dτ c (τ )

(u · ∇ψ)ρ dv +

c (τ )

(u · n f )ρ ψ ds.

∂ c (τ )

A semi-discrete formulation is obtained by replacing functions ρ and ψ addition, to enforce the conservation law, we introduce a unique value definition of a mass flux between cells. Finally, since velocity u is well mesh cells, we approximate it there using specially designed projectors



d





ρh ψh dv = −

dτ c (τ )

 c,k+1 (u) · ∇ψh ρh dv +

c (τ )



by their cell-based Taylor expansions ρh and ψh . In of function ρ on face f , which leads to the unique defined only on mesh faces and is unknown inside c,k+1 (u), k ≥ 0, described later:

(u · n f )ρh∗ ψh ds.

f (τ )

Now, we sum up these cell-based equations and change the summation order for the interface terms from cells to faces. Since u · n f = 0 on the domain boundary, the last summation below uses only internal faces:

 d 

dτ

c

ρh ψh dv = −

  c

c (τ )

  c,k+1 (u) · ∇ψh ρh dv +



(u · n f )ρh∗ [ψh ] ds.

(4)

f f (τ )

c (τ )

Since ψh is a polynomial only at pseudo-time τ = 0, we have to perform integration in the reference coordinate system ξ . Let matrix J denote the Jacobian of the space-time map x:

J ik =

∂ xi ∂ Fi =τ + (1 − τ )δik , ∂ξk ∂ξk

j = det(J),

where i , k = 1, 2. We recall the Nanson formula for the transformation of normal vectors and the change of variables formula for the gradient operator:

n f ds = j J− T N f dS ,

∇x ρ = J−T ∇ξ ρ ,

where dS and ds are surface measures in ξ and x coordinate systems. Since, we do not have an explicit formula for J inside mesh cells, we approximate the Jacobian using the available approximation of velocity u. Formula (1) leads to





Jh = I + τ ∇ξ c ,k+1 (u) ,

j h = det(Jh ),

where I is the identity matrix. Inserting the above formulas in (4), we obtain the equivalent weak formulation in the reference coordinate system:

 d 

dτ

c

ρh ψh jh dV = −

  c

c ( 0)

 c,k+1 (u) · jh Jh−T ∇ξ ψh ρh dV

c ( 0)

+



(u · j J−T N f )ρh∗ [ψh ] dS .

(5)

f f ( 0)

The stability argument suggests to select the downwind value for

1

(u · n f )ρh∗ = (u · n f ){ρh } − |u · n f | [ρh ]. 2

ρh∗ , e.g. (6)

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

5

For high-order DG schemes, flux u · n f may change sign on face f . Numerical experiments show that in the extreme case of a random motion of mesh vertices, this happens on 40% of mesh edges. In practice, the face integral is calculated using a quadrature rule and formula (6) provides flux values at the quadrature points. The integrability argument suggests to define c ,k+1 as a projector onto a polynomial space of order k + 1, see Section 4.1 for more detail. Recall that j J− T is called in linear algebra the matrix of cofactors for J. Since ∇ξ c ,k+1 (u) is a matrix of polynomials, the entries of the matrix of cofactors for Jh are also polynomials, and j h is a polynomial as well. In such a case, the divergence theorem can be used to reduce the volume integral to face integrals for any polygonal cell c. For an efficient practical implementation, we recommend to combine the divergence theorem with Euler’s homogeneous function theorem [26,14]. The vector function j J− T N f , related to the face normal, is defined completely by the current position of face f , which is important for calculating a unique flux across this face. Recall that the map F(ξ ) for each face is a vector polynomial by our assumption, and so is the remap velocity u. Let θ be the natural parameterization of curved face f . Then, u · j J− T N f dS = ϕ (θ) dθ , where ϕ (θ) is a polynomial. Thus, all integrals in (5) are computable. 4. Fully discrete scheme We use the third order TVD Runge-Kutta (RK) scheme [19] to derive a fully discrete scheme. First, we rewrite the weak formulation in a matrix form:

d dt

(Mρ h ) = G(ρ h ),

where ρ h is the algebraic vector of coefficients in all Taylor expansions, M represents the volumetric matrix (the mass matrix in a finite-element language) and G is the discrete functional combining volumetric and surface terms. Second, we change variable φ h := Mρ h to obtain a conventional system of ODE equations. 4.1. Virtual element projector on a polygonal cell To complete the formulation of the discrete scheme, we need to define projector c ,k+1 . We use superscript k + 1 (instead of k) to stress the fact that we build a scheme of order k + 1. In other words, given a cell c and a velocity function on mesh faces, we want to build its polynomial approximation inside c. The analysis presented in [30] shows that projection of velocity u on a polynomial space made with respect to the energy semi-norm has some benefits, see also Remark 4.1 below. We start with building the auxiliary projector c0,k+1 using the least-square algorithm. We assume that this projector is well defined by the velocity on cell boundary. In this case, it becomes the L 2 projector for the following virtual element space:

 Vk+1 (c ) = v ∈ C 0 (¯c ) : v| f ∈ P k ( f ) ∀ f ∈ ∂ c , v ∈ P k+1 (c ),    v q dV = c0,k+1 (v) q dV ∀q ∈ P k+1 (c ) . c

c

Hereafter, P denotes the space of polynomials (or vector polynomials) of order at most l. The degrees of freedom for this space are velocity values at cell vertices and velocity moments of order k − 2 on cell faces. Indeed, consider function u ∈ Vk+1 (c ) and assume that all degrees of freedom for it are zeros. The definition of the virtual space states that u is a polynomial of order at most k on each edge f . Now, k − 2 zero moments and two zero values at the edge end points imply that u = 0 on f . By our assumption, the L 2 projector is defined uniquely by the boundary degrees of freedom. Thus c0,k+1 (u) = 0. Finally, using this and u ∈ P k+1 (c ) in the integration by parts formula, we obtain l





|∇ u|2 dV = − c



u · u dV + c



(∇ u) · u dS = − ∂c



u · c0,k+1 (u) dV = 0.

u · u dV = − c

c

Hence, u is a constant vector function which is zero on ∂ c. We conclude that u = 0 which shows uniqueness. The existence of a virtual function and the unisolvency property follow from the theory of elliptic operators, see [30] for more detail. We define projector c ,k+1 (u) : Vk+1 (c ) → P k+1 (c ) as the solution of





∇ c,k+1 (u) : ∇ q dV = c

c

subject to





c,k+1 (u) dV = c

∇ u : ∇ q dV

u dV . c

∀q ∈ P k+1 (c ),

(7)

6

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

Equation (7) defines a polynomial up to a constant, which is fixed using the constraint. Integrating by parts the right-hand side of (7) and using the definition of the virtual space, we obtain the equivalent definition of the projector:





c



c0,k+1 (u) · q dV +

∇ c,k+1 (u) : ∇ q dV = − c

u · (∇ q · n) dS

∀q ∈ Pk+1 (c ).

(8)

∂c

Note that all integrals are computable. Remark 4.1. The virtual element space is continuous across mesh faces; hence, it is more suitable for analysis, while the projector is more suitable for a computer implementation. Taking q = c ,k+1 (u) in (7), we obtain:



∇ c,k+1 (u) 2L 2 (c) =

∇ u : ∇ c,k+1 (u) dV ≤ ∇ c,k+1 (u) L 2 (c) ∇ u L 2 (c) ,

(9)

c

i.e. this projector reduces the energy norm of the resulting polynomial. Thus, this polynomial approximation of u does not have spurious oscillations. Note that the auxiliary projector does not have this property. Remark 4.2. If c is a square, then c ,k+1 (u) preserves bilinear and biquadratic functions. This establishes a connection with the theory of the finite element maps. 4.2. The conservation law It is easy to see that the choice ψh = 1 in formula (5) leads to the conservation of the weighted integral of

 d 

dτ

c

ρh jh dV = 0.

c ( 0)

Using the orthogonality property of basis functions at time this equation as follows:

 c

ρh :

ρhn+1 jh dV =

τ = 0 and the fact that jh = 1 at this time moment, we rewrite

 (ρ0c )n |c (0)|. c

c ( 0)

Note the presence of jh in the left-hand side. Since, the Lagrangian step of a simulation relies on the fact that mean values are responsible for the conservation law, a re-initialization procedure is needed before starting the next Lagrangian step. The projection procedure described below takes j h into account. For each cell c (0) ∈ h , the remap velocity u defines explicitly trajectories of boundary points. Since each basis function ˜ h and is a remains constant along these trajectories, the remapped solution is known on the boundary of cell c (1) ∈  polynomial on its face. Thus, we can project the non-polynomial remapped solution, ρhN , on the space of polynomials using the scalar version of the energy projector, e.g.



˜ c,k (ρ N ) · ∇ q dV = ∇ h

c (1 )

∀q ∈ P k (c (1)),

∇ ρh · ∇ q dV

(10)

c (1 )

subject to





˜ c,k (ρ N ) dV =  h

c (1 )



ρhN jh dV .

c ( 0)

The last formula ensures the conservation law. 4.3. Limiting of discontinuous solutions Remaps of a discontinuous solution requires a limiting strategy to reduce significantly oscillations around discontinuities. As the proof of principle, we use the Barth-Jespersen limiter, although we see no obvious obstacles in using accuracy preserving limiters such as the hierarchical limiter from [25]. Let ρˆh denote the limited solution. Specifics of the remap problem is that the limiting strategy consists of two steps. In the first step, we limit solution ρh in each cell by multiplying all high-order terms by a single non-negative factor:







ρˆh c = ρ0c + αc ρ1c ψ1c + ρ2c ψ2c + ρ3c ψ3c + . . . , 0 ≤ αc ≤ 1.

(11)

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

7

Fig. 2. Top row: original (left) and two modified (right) quadrilateral meshes. Bottom row: original (left) and two modified (right) polygonal meshes.

αc is calculated such that the limited function satisfies the following inequalities: 



k min {ρ0c } ≤ ρˆh c (ξ k ) ≤ max {ρ0c }, ∀ξ k ∈ Gc 1 ,



The factor

c ∈Fc

c ∈Fc

(12)

k

where Fc is the set of cells that have at least one common point with cell c and Gc 1 is a set of limiting points. The value k1 = k + 1 that we use in numerical experiments corresponds to the Gauss-Legendre quadrature rule of degree 2k + 1. We think that this quadrature rule is sufficiently accurate to be used in the future for the scheme optimization. Analysis of an optimal quadrature rule is beyond the scope of this work. The conservation law implies that the integral of ρh jh must be preserved during the limiting strategy. Simple calculations show that this is achieved by shifting the mean value of the limited function in each mesh cell. Thus, the second step of the limiting strategy is

ρˆ0c = αc ρ0c + (1 − αc )

1

|c (τ )|



ρh jh dV .

(13)

c (τ )

5. Numerical experiments

˜ h , we consider two mesh relaxation processes defined by the Taylor-Green vorticial To generate a modified mesh  motion (TG map) and the compression/expansion map (CE map). In the TG map, the nodes of a given mesh are moved according to the following ODEs: d xˆ dτ

= 0.2 sin(π x) cos(π y ),

d yˆ dτ

= −0.2 cos(π x) sin(π y ),

with the initial values xˆ (0) := ξ and the final dimensionless time of 1. Note that the remap velocity u linearizes this motion and thus differs from dxˆ /dτ . The CE map is given by the cubic velocity:

ˆ= u

1 2

ξ1 ξ2 (1 − ξ ).

We consider two sequences of quadrilateral and general polygonal meshes. The original and two modified meshes are shown in Fig. 2 for modest mesh resolutions. Calculation of the discrete L 2 error, denoted by  (ρh ), deserves a few comments. Since the remap velocity u is a virtual function, in general, it has no well-defined point values inside a mesh cell. For this reason, we cannot calculate a trajectory for any interior point. In contrast, trajectories of boundary points are known. For instance, we can use values of basis functions at vertices ξ v for the error calculation:

 (ρh ) =

 |c | 

1/2 |ρ (˜x v ) − ρh (˜x v ))|2 , Nc v ∈c ˜h c ∈

where x˜ v is the position of vertex v on the modified mesh, and Nc is the number of vertices in cell c.

(14)

8

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

Fig. 3. TG map on square meshes: solution isolines in three DG(P k ) schemes with k = 0 (left panel), k = 1 (middle panel) and k = 2 (right panel).

˜ c,k (ρ N ) which is a polynomial in cell Remark 5.1. Another way to estimate the true error is to use the projected solution  h c. Our numerical experiments with smooth solutions indicate that the projection step changes the discrete L 2 error (14) by about 1% and the maximum error by 3% to 8%. The visualization is done with the VisIt software [13] using the solution mean values ρ0c . Although we use curved mesh edges in a computer program for k > 1, they are visualized using straight lines in the presented figures. 5.1. Numerical integration over polygonal cells Efficiency of the proposed technology is based on the feasibility to integrate polynomial functions over arbitrarily-shared cells. For this purpose, we recall the Euler homogeneous function theorem. Let f (ξ ) be a continuously differential homogeneous function of degree α :

f (λξ ) = λα f (ξ ),

α > 0.

This function satisfies the Euler homogeneous function theorem:

α f (ξ ) = ξ · ∇ξ f (ξ ). The Gauss-Green formula holds for any polygonal cell c:





divξ ξ f (ξ ) dV + c

c

Note that divξ ξ = 2. Thus,



f (ξ ) dV = c



ξ · ∇ξ f (ξ ) dV =

1 2+α

(ξ · n) f (ξ ) dS . ∂c



(ξ · n) f (ξ ) dS .

f ∈∂ c f

Every polynomial q(ξ ) is a linear combination of monomials that are homogeneous functions. Since ξ · n is a polynomial with respect to the edge parameterization, the integral of q(ξ ) over cell c is reduced to the problem of integrating polynomials over its faces f . The later could be performed efficiently using the Gauss-Legendre quadrature formulas. 5.2. Smooth solution: logically square meshes Let  be the unit square. To study the formal order of convergence of the proposed scheme, we consider a smooth solution given by

ρ (ξ1 , ξ2 ) = sin(6ξ1 ) sin(3ξ2 ).

(15)

In the first experiment, we compare the unlimited DG(P ) solutions for the TG map. The results are collected in Table 1 and in Fig. 3. Hereafter, the convergence rates are calculated using the linear regression algorithm. We suppress a time integration error by using the time step that is smaller than that required for scheme stability. The time step on the coarsest mesh is 0.01 and is reduced twice for each mesh refinement. Notice that with each increment of k, we gain almost one order of accuracy on the coarsest mesh and almost two orders on the finest mesh. The DG(P 2 ) scheme uses quadratic approximation of edges for τ > 0. In some experiments, the optimal order of convergence holds also for linear meshes. In general, the lack of an accurate approximation of the remap velocity on mesh edges k

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

9

Table 1 TG map: the discrete L 2 and maximum errors for DG(P k ) schemes on a sequence of square meshes. 1/h

DG(P 0 )

DG(P 1 )

DG(P 2 )

L2

L∞

L2

L∞

L2

L∞

16 32 64 128

1.125e-1 5.772e-2 2.924e-2 1.471e-2

2.279e-1 1.317e-1 6.699e-2 3.387e-2

1.610e-2 4.095e-3 1.032e-3 2.592e-4

4.492e-2 1.182e-2 3.032e-3 7.677e-4

1.093e-3 1.346e-4 1.670e-5 2.085e-6

3.553e-3 4.828e-4 6.305e-5 8.047e-6

rate

0.979

0.922

1.986

1.957

3.011

2.930

Table 2 CE map: the discrete L 2 and maximum errors for DG(P k ) schemes on a sequence of square meshes. 1/h

DG(P 0 ) L

DG(P 1 )

2

L∞

L

2

DG(P 2 ) L∞

L2

L∞

16 32 64 128

1.123e-1 5.715e-2 2.881e-2 1.446e-2

2.678e-1 1.416e-1 7.215e-2 3.634e-2

1.650e-2 4.196e-3 1.056e-3 2.646e-4

4.404e-2 1.169e-2 2.988e-3 7.552e-4

1.149e-3 1.395e-4 1.733e-5 2.171e-6

3.699e-3 5.049e-4 6.561e-5 8.351e-6

rate

0.986

0.962

1.988

1.957

3.015

2.932

Table 3 TG map: the discrete L 2 and maximum errors for DG(P k ) schemes on a sequence of polygonal meshes. #cells

DG(P 0 ) L

2

2

DG(P 1 ) L∞

L

2

DG(P 2 ) L∞

L2

L∞

16 322 642 1282

1.430e-1 8.134e-2 4.671e-2 2.676e-2

4.305e-1 3.572e-1 2.725e-1 1.993e-1

1.977e-2 4.694e-3 1.139e-3 2.805e-4

7.848e-2 1.944e-2 4.662e-3 1.165e-3

1.568e-3 1.791e-4 2.128e-5 2.608e-6

9.311e-3 1.150e-3 1.378e-4 1.842e-5

rate

0.805

0.372

2.046

2.028

3.077

3.001

Table 4 CE map: the discrete L 2 and maximum errors for DG(P k ) schemes on a sequence of polygonal meshes. #cells

DG(P 0 )

DG(P 1 )

DG(P 2 )

L2

L∞

L2

L∞

L2

L∞

162 322 642 1282

1.230e-1 6.434e-2 3.352e-2 1.753e-2

4.153e-1 2.253e-1 1.644e-1 1.233e-1

2.026e-2 4.920e-3 1.209e-3 3.007e-4

8.576e-2 2.865e-2 7.241e-3 1.814e-3

1.666e-3 1.944e-4 2.413e-5 3.063e-6

1.401e-2 1.819e-3 2.207e-4 2.774e-5

rate

0.937

0.571

2.025

1.867

3.027

2.998

leads to the loss of the optimal convergence rate. A similar observation was made in [30] regarding DG schemes for a level set problem. In the second experiment, we compare the DG(P k ) solutions for the CE map. The results collected in Table 2 confirm our previous conclusions. In the case of square meshes, the projector c ,k+1 recovers the finite element map. Our results are in agreement with the results reported in [3] for the finite element method on curvilinear meshes. 5.3. Smooth solution: polygonal meshes In this section, we consider again the smooth solution (15) but the sequence of polygonal meshes shown in Fig. 2. In two sets of experiments, we compare the DG(P k ) solutions for the TG and CE maps. The results are collected in Tables 3 and 4. The remapped solutions are shown in Fig. 4. In contrast to the square meshes, the projected velocity c ,k+1 (u) is discontinuous across polygonal cells. Recall, that the virtual element space is continuous by construction. The benefits of using higher-order DG schemes are evident also in the case of polygonal meshes. The lack of optimal convergence of the DG(P 0 ) scheme requires a rigorous analysis. We think that it could be connected to the geometric conservation law which was used implicitly in the derivation of the conservative weak formulation. Let us compare the DG schemes for a given level of accuracy. On a powerful OSX computational platform, the DG(P 0 ) scheme on the finest mesh runs approximately 200 times longer than the DG(P 1 ) scheme of the coarsest mesh and still

10

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

Fig. 4. TG map on polygonal meshes: solution isolines in three DG(P k ) schemes with k = 0 (left panel), k = 1 (middle panel) and k = 2 (right panel).

Fig. 5. The initial states on unstructured quadrilateral meshes with mesh resolutions h = 1/40 (left panel) and h = 1/80 (right panel).

gives a larger error. The DG(P 1 ) scheme on the mesh with mesh resolution 1/64 runs approximately 14 times longer than the DG(P 2 ) scheme of the coarsest mesh. 5.4. Discontinuous solution Shock structures are typical for hyperbolic systems. In this section, we consider the remap of a solution that combines discontinuous and continuous parts. The problem setup is based on the solid body rotation test described in [28]. We also study the impact of the Barth-Jespersen limiter on the smooth and non-smooth parts of this solution. The computational domain is the disk of unit radius centered at the origin. The discontinuous part of the solution is given by the characteristic function of the notched disk [36,45]. The disk has radius r0 = 0.4 and is centered at point (0.5, 0). The notch is formed by removing the intersection of the disk with the rectangle of width r0 /4 and height r0 , see Fig. 5. The continuous part of the solution with discontinuous first derivatives is given by the cone of radius r1 = 0.35 centered at point (cos(φ)/2, sin(φ)/2), where φ = 120◦ . The cone has the unit height. The smooth part of the solution is given by the hump that has circular support of radius r2 = 0.35 centered at point ξ 2 = (cos(φ)/2, − sin(φ)/2). The hump shape is given by

ρ (ξ ) =

1 4

(1 + cos(π ξ − ξ 2 /r2 )).

We consider two unstructured quadrilateral and two polygonal meshes with effective mesh resolutions 1/40 and 1/80, see Figs. 5 and 8. The color map corresponds to the solution value. The solution is non-negative. Thus, to visualize clearly the underlying mesh, we color the solution only when its value is above 0.01. We rotate each mesh h counter clockwise 360◦ in increments of 9◦ . On each segment of rotation, we calculate the remap velocity and remap solution using the method described above. In total, we perform 40 advection solves. The time step on the coarse and fine meshes are 0.025 are 0.0125, respectively. In Figs. 6 and 7 we compare limited and unlimited solutions on two unstructured quadrilateral meshes. In all experiments, the conservation law is satisfied up to machine precision. The lack of limiters leads to moderate overshoots and undershoots around the solution discontinuity. The violation of bounds is between 0.05 and 0.11, depending on the experiment, and is slightly less for the DG(P 2 ) solution. The Barth-

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

11

Fig. 6. The DG(P 1 ) (left column) and DG(P 2 ) (right column) solutions after 360◦ rotation of the unstructured quadrilateral mesh with the mesh resolution h = 1/40. Top row shows the unlimited solutions. Bottom row shows the limited solutions.

Fig. 7. The DG(P 1 ) (left column) and DG(P 2 ) (right column) solutions after 360◦ rotation of the quadrilateral mesh with the mesh resolution h = 1/80. Top row shows the unlimited solutions. Bottom row shows the limited solutions.

12

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

Fig. 8. The initial states on polygonal meshes with mesh resolutions h = 1/40 (left panel) and h = 1/80 (right panel).

Fig. 9. The DG(P 1 ) (left column) and DG(P 2 ) (right column) solutions after 360◦ rotation of the polygonal mesh with the mesh resolution h = 1/40. Top row shows the unlimited solutions. Bottom row shows the limited solutions.

Jespersen limiter removes the oscillations. The reported violation of bounds is between 0 and 0.0003, due to the fact that the scheme is formally not maximum principle preserving. As expected, the Barth-Jespersen limiter smears the cone pike. A possible remedy is to use shock detectors, e.g. like that in [32] together with accuracy-preserving limiters, e.g. like that in [25]. Our preliminary work revealed that the shock detectors developed for the Lagrangian hydrodynamics require additional analysis for the remap problem which is beyond the scope of this work. The DG(P 2 ) scheme seems to be less diffusive and preserves better the shape of the notched cylinder than the DG(P 1 ) scheme. It also reduces smearing of the cone tip. At the same time, it leads to more rounded brims of the notched cylinder. This indicates that higher-order limiters may be needed for higher-order DG schemes. In Figs. 9 and 10 we compare limited and unlimited solutions on two unstructured polygonal meshes. The conservation law is satisfied again up to machine precision. Our previous conclusions regarding overshoots in the unlimited solutions and the impact of the limiter hold in this case as well.

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

13

Fig. 10. The DG(P 1 ) (left column) and DG(P 2 ) (right column) solutions after 360◦ rotation of the polygonal mesh with the mesh resolution h = 1/80. Top row shows the unlimited solutions. Bottom row shows the limited solutions.

6. Conclusion The remap problem is a part of many arbitrary Lagrangian-Eulerian (ALE) hydrodynamics simulations. We proposed a new technology to remap discontinuous Galerkin solutions between two polygonal meshes with curved edges. The technology is based on solving a linear advection equation and is extendable to three dimensions. A map between two meshes is known only on cell boundaries and is extended inside cells using the recently introduced virtual element method. This method produces a remap velocity that lives in a continuous virtual element space. In the proposed scheme, this uncomputable velocity is projected on a polynomial space in each mesh cell using the energy-based virtual projectors. For smooth solutions, we demonstrated the k + 1 order of convergence of the DG(P k ), k = 1, 2, scheme on square and polygonal meshes in both discrete L 2 and maximum norms. The numerical experiments have shown significant computational advantages for using these schemes compared to the DG(P 0 ) scheme. For discontinuous solutions, we proposed a new conservative limiter based on the scalar Barth-Jespersen limiter. The numerical experiments on quadrilateral and polygonal meshes have shown significant reduction of oscillations in the limited solution. 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. Acknowledgements This work was carried out under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396 and the Laboratory Directed Research and Development (LDRD) Program. We thank Dr. Xiaodong Liu for useful discussions on limiters and shock detectors for DG schemes. We thank Mack Kenamond for his help in using computer software that builds polygonal meshes. References [1] V. Aizinger, C. Dawson, B. Cockburn, P. Castillo, Local discontinuous Galerkin method for contaminant transport, Adv. Water Resour. 24 (2000) 73–87.

14

K. Lipnikov, N. Morgan / Journal of Computational Physics 399 (2019) 108931

[2] R. Anderson, V. Dobrev, T. Kolev, D. Kuzmin, M. Quezada de Lun, R. Rieben, V. Tomov, High-order local maximum principle preserving (MPP) discontinuous Galerkin finite element method for transport equation, J. Comput. Phys. 334 (2017) 102–124. [3] R. Anderson, V. Dobrev, T. Kolev, R. Rieben, Monotonicity in high-order curvilinear finite element arbitrary Lagrangian-Eulerian remap, Int. J. Numer. Methods Fluids 77 (2015) 249–273. [4] B. Ayuso de Dios, K. Lipnikov, G. Manzini, The nonconforming virtual element method, M2AN: Math. Model. Numer. Anal. 50 (3) (2016) 879–904. [5] F. Baaijens, Application of low-order discontinuous Galerkin methods to the analysis of viscoelastic flows, J. Non-Newton. Fluid Mech. 52 (2) (1994) 37–57. [6] T. Barth, D. Jespersen, The design and application of upwind schemes on unstructured meshes, in: 27th Aerospace Sciences Meeting, AIAA-1989-0366, Reno, NV, 1989. [7] F. Bassi, S. Rebay, High-order accurate discontinuous finite element solution of the 2D Euler equations, J. Comput. Phys. 138 (1997) 251–285. [8] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L.D. Marini, A. Russo, Basic principles of virtual element methods, Math. Models Methods Appl. Sci. 23 (2013) 119–214. [9] L. Beirao da Veiga, K. Lipnikov, G. Manzini, The Mimetic Finite Difference Method for Elliptic PDEs, Springer, 2014, 408 pp. [10] D. Benson, An efficient, accurate, simple ALE method for nonlinear finite element programs, Comput. Methods Appl. Mech. Eng. 72 (1989) 305–350. [11] E. Caramana, D. Burton, M. Shashkov, P. Whalen, The construction of compatible hydrodynamics algorithms utilizing conservation of total energy, J. Comput. Phys. 146 (1998) 227–262. [12] Z. Chen, B. Cockburn, C. Gardner, J. Jerome, Quantum hydrodynamic simulation of hysteresis in the resonant tunneling diode, J. Comput. Phys. 117 (1995) 274–280. [13] H. Childs, E. Brugger, B. Whitlock, J. Meredith, S. Ahern, D. Pugmire, K. Biagas, M. Miller, C. Harrison, G. Weber, H. Krishnan, T. Fogal, A. Sanderson, C. Garth, E. Wes Bethel, D. Camp, O. Rübel, M. Durant, J. Favre, P. Navrátil, VisIt: an end-user tool for visualizing and analyzing very large data, in: High Performance Visualization–Enabling Extreme-Scale Scientific Insight, Oct 2012, pp. 357–372. [14] E. Chin, J. Lasserre, N. Sukumar, Numerical integration of homogeneous functions on convex and nonconvex polygons and polyhedra, Comput. Mech. 56 (6) (2015) 967–981. [15] B. Cockburn, C.-W. Shu, The Runge-Kutta discontinuous Galerkin finite element method for conservation laws V: multidimensional systems, J. Comput. Phys. 141 (1998) 199–224. [16] V. Dobrev, T. Kolev, R. Rieben, High-order curvilinear finite element methods for Lagrangian hydrodynamics, SIAM J. Sci. Comput. 34 (5) (2012) 606–641. [17] J. Dukowicz, J. Baumgardner, Incremental remapping as a transport/advection algorithm, J. Comput. Phys. 160 (2000) 318–335. [18] J. Dukowicz, J. Kodis, Accurate conservative remapping (rezoning) for arbitrary Lagrangian-Eulerian computations, SIAM J. Stat. Comput. 8 (1987) 305–321. [19] S. Gottlieb, C.-W. Shu, Total variation diminishing Runge-Kutta schemes, Math. Comput. 67 (221) (1998) 73–85. [20] P. Gremaud, J. Matthews, Discontinuous Galerkin methods for Friedrichs’ systems. I. General theory, J. Comput. Phys. 166 (2) (2001) 63–83. [21] C. Hirt, A. Amsden, J. Cook, An arbitrary Lagrangian-Eulerian computing method for all flow speeds, J. Comput. Phys. 14 (1974) 227–253. [22] Z. Jibben, M. Herrmann, An arbitrary-order Runge-Kutta discontinuous Galerkin approach to reinitialization for banded conservative level sets, J. Comput. Phys. 349 (2017) 453–473. [23] D. Kershaw, M. Prasad, M. Shaw, J. Milovich, 3D unstructured mesh ALE hydrodynamics with the upwind discontinuous finite element method, Comput. Methods Appl. Mech. Eng. 158 (1998) 81–116. [24] M. Kucharık, M. Shashkov, B. Wendroff, An efficient linearity-and-bound-preserving remapping methods, J. Comput. Phys. 188 (2003) 462–471. [25] D. Kuzmin, A vertex-based hierarchical slope limiter for p-adaptive discontinuous Galerkin methods, J. Comput. Phys. 233 (12) (2010) 3077–3085. [26] J. Lasserre, Integration and homogeneous functions, Proc. Am. Math. Soc. 127 (3) (1999) 813–818. [27] O. Lepsky, C. Hu, C.-W. Shu, Analysis of the discontinuous Galerkin method for Hamilton-Jacobi equations, Appl. Numer. Math. 33 (2000) 423–434. [28] R. LeVeque, High-resolution conservative algorithms for advection in incompressible flow, SIAM J. Numer. Anal. 33 (1996) 627–665. [29] Z. Li, X. Yu, Z. Jia, The cell-centered discontinuous Galerkin method for Lagrangian compressible Euler equations in two-dimensions, Comput. Fluids 96 (2014) 152–164. [30] K. Lipnikov, N. Morgan, A High-Order Discontinuous Galerkin Method for Level Set Problems on Polygonal Meshes, Technical Report 2018-18-29017, Los Alamos National Laboratory, 2018. [31] X. Liu, N. Morgan, D. Burton, A Lagrangian discontinuous Galerkin hydrodynamic method, Comput. Fluids 163 (2018) 68–85. [32] X. Liu, N. Morgan, D. Burton, A high-order Lagrangian discontinuous Galerkin hydrodynamic method for quadratic cells using a subcell mesh stabilization scheme, J. Comput. Phys. 386 (1) (2019) 110–157. [33] R. Loubére, P.-H. Maire, M. Shashkov, J. Breil, S. Galera, ReALE: a reconnection-based arbitrary Lagrangian-Eulerian method, J. Comput. Phys. 229 (2010) 4724–4761. [34] R. Loubére, J. Ovadia, R. Agrall, A Lagrangian discontinuous Galerkin-type method on unstructured meshes to solve hydrodynamics problems, Int. J. Numer. Methods Fluids 44 (2004) 645–663. [35] L. Margolin, M. Shashkov, Second-order sign-preserving conservative interpolation (remapping) on general grids, J. Comput. Phys. 184 (1) (2003) 266–298. [36] N. Morgan, J. Waltz, 3D level set methods for evolving fronts on tetrahedral meshes with adaptive mesh refinement, J. Comput. Phys. 336 (2017) 492–512. [37] N. Nguyen, J. Peraire, B. Cockburn, An implicit high-order hybridizable discontinuous Galerkin method for the incompressible Navier-Stokes equations, J. Comput. Phys. 230 (4) (2011) 1147–1170. [38] J. Peery, D. Carroll, Multi-material ALE methods in unstructured grids, Comput. Methods Appl. Mech. Eng. 187 (2000) 591–619. [39] W. Reed, T. Hill, Triangular Mesh Methods for the Neutron Transport Equation, Technical Report LA-UR-73-479, Los Alamos National Laboratory, 1973. [40] C.-W. Shu, J. Yan, A local discontinuous Galerkin method for KdV type equations, SIAM J. Numer. Anal. 40 (2002) 769–791. [41] B. Taylor, Methodus incrementorum directa et inversa, Impensis Gulielmi Innys, Londini, 1717, 119 pp. [42] W.F. Trench, B. Kolman, Multivariable Calculus with Linear Algebra and Series, Academic Press, 1972, 758 pp. [43] F. Vilar, Cell-centered discontinuous Galerkin discretization for two-dimensional Lagrangian hydrodynamics, Comput. Fluids 64 (2012) 64–73. [44] F. Vilar, P.-A. Maire, R. Abgrall, A discontinuous Galerkin discretization for solving two-dimensional gas dynamics equations written under total Lagrangian formulation on general unstructured grids, J. Comput. Phys. 276 (2014) 188–234. [45] S. Zalesak, Fully multidimensional flux-corrected transport algorithms for fluids, J. Comput. Phys. 31 (1979) 335–362. [46] X. Zhang, C.-W. Shu, On positivity-preserving high-order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes, J. Comput. Phys. 229 (2010) 8918–8934.