Curl constraint-preserving reconstruction and the guidance it gives for mimetic scheme design

Several important PDE systems, like magnetohydrodynamics and computational electrodynamics, are known to support involutions where the divergence of a vector field evolves in divergence-free or divergence constraint-preserving fashion. Recently, new classes of PDE systems have emerged for hyperelasticity, compressible multiphase flows, so-called first order reductions of the Einstein field equations, or a novel first order hyperbolic reformulation of Schr\"odinger's equation, to name a few, where the involution in the PDE supports curl-free or curl constraint-preserving evolution of a vector field. Since mimetic numerical schemes for the solution of the former class of PDEs are well-developed, we draw guidance from them for the solution of the latter class of PDEs. We show that a study of the curl constraint-preserving reconstruction gives us a great deal of insight into the design of consistent, mimetic schemes for these involutionary PDEs. The importance of multidimensional Riemann solvers in facilitating the design of such schemes is also documented. We study the problem of curl constraint-preserving reconstruction as it pertains to the design of mimetic discontinuous Galerkin (DG) and finite volume (FV) schemes for PDEs that support such an involution. This is done for two and three dimensional structured mesh problems where we deliver closed form expressions for the reconstruction. The role that this reconstruction plays in the curl-free, or curl-preserving prolongation of vector fields in adaptive mesh refinement (AMR) is also discussed. In two dimensions, a von Neumann analysis of structure-preserving DG-like schemes that mimetically satisfy the curl constraints, is also presented. Numerical results are also presented to show that the schemes meet their design accuracy.


I.1) Introduction
There has been a lot of emerging interest in mimetic scheme design. These are schemes that preserve structures in the solution that arise from involutions in the governing PDEs. In other words, the PDE itself has some extra symmetries that result in certain features of the solution remaining invariant; and we want the numerical scheme to mimic that.
The simplest example of such involution-constrained PDEs consists of the magnetohydrodynamic (MHD) equations, where Faraday's law ensures divergence-free evolution of the magnetic induction vector field. Another prominent example consists of computational electrodynamics (CED) -the numerical solution of Maxwell's equations -where the divergences of the magnetic induction vector field and the electric displacement vector field are held zero as long as radiation does not interact with a conductor. Numerous papers have been written on these topics, where it has been realized that the divergence-preserving reconstruction of vector fields is an important building block for scheme design and adaptive mesh refinement (AMR) (Balsara and Spicer [3], Balsara [5], [6], [8], Balsara and Dumbser [14], Xu et al. [61], Balsara and Käppeli [20], [24], Balsara et al. [10], [18], [22], [23], [25], [26] and Hazra et al. [44]). To get fully constraint-preserving, mimetic, time-evolution it was realized that certain update variables have to be collocated at certain favored locations on a mesh. To obtain such variables in a properly upwinded fashion, it is crucially important to invoke a multidimensional Riemann solver at the edges of the computational mesh (Balsara, [9], [10], [13], Balsara, Dumbser and Abgrall [12], Balsara and Dumbser [15], Balsara et al. [17], Balsara and Nkonga [21]). The multidimensional Riemann solver is, therefore, the other building block of the scheme. In obtaining highly accurate globally constraint-preserving discontinuous Galerkin-like (DG-like) schemes for MHD and CED, Balsara and Käppeli [20], [24] showed that both these building blocks were crucially important. They showed that if one attempts to bypass either of these building blocks, it will result in an unstable DG-like scheme. We call these schemes DG-like because they evolve all the facebased modes of the vector field so as to ensure divergence constraint-preserving evolution of the vector field; however, they are not like classical DG schemes because the modes are not defined on the volumes. We, therefore, see that a study of the involution-preserving reconstruction can provide substantial insights into scheme design. A review of globally divergence constraintpreserving DG schemes for CED is also available in Balsara and Simpson [27] which collects all the ideas together in an easily accessible format in one place.
While MHD and CED are relatively well-studied PDEs with a divergence constraint, a new class of PDEs has recently emerged and their involution constraints are equally interesting. We are referring to PDEs that support curl-free evolution of vector fields. Indeed, the evolution is curlfree in these systems only as long the source terms in the governing equations are zero. Numerous PDEs of great practical interest fall in this category. Many of the hyperbolic systems resulting from the Godunov-Peshkov-Romenski (GPR) formulation for hyperelasticity and compressible multiphase flow with and without surface tension have such curl-preserving update equations (Godunov and Romenski [42], Romenski [51], Romenski et al. [52], Peshkov and Romenski [49], Dumbser et al. [36], [37], Schmidmayer et al. [53]). The equations of General Relativity when cast in the FO-CCZ4 formulation also have such a structure (Alic et al. [1], [2], 2012, Brown et al. [32], Dumbser et al. [38], Dumbser, et al. [39]). Similarly, it has recently become possible to recast Schrödinger's equation in first order hyperbolic form, and the time-evolution of this very important equation also has curl-preserving constraints (Dhaouadi et al. [35]). As with the divergencepreserving reconstruction, the curl-preserving reconstruction also plays an important role in guiding scheme design. The goal of this paper is to show how curl-preserving reconstruction of vector fields can be carried out and why it is so important in the design of curl-constraintpreserving schemes. We restrict our focus to structured meshes, since the treatment of unstructured meshes will be the topic of another paper.
In this work we take on the task of designing a globally curl constraint-preserving reconstruction. This means that the curl of a vector field, evaluated over any closed loop, is always either zero or equal to a specified divergence-free vector field. It may even prove advantageous to

I.2) Introduction : Intercomparing Divergence-preserving and Curl-Preserving Schemes
Let us take a look at a representative divergence-preserving system (we pick Maxwell's equations) and compare and contrast it with the simplest curl-preserving system (we pick the equations for a fluid with hyperbolic heat conduction). We hope to show that there are important points of similarity and also some differences. This compare and contrast will be instructive in terms of the kinds of constraint-preserving schemes that we design for PDE systems with involution constraints.
For Maxwell's equations, the primal variables are the electric displacement D and the magnetic induction B which evolve according to the curl of the magnetic field H and the curl of the electric field E . The constitutive relations between these vector fields are = D εE and = B μH where ε is the permittivity tensor and μ is the permeability tensor. Because of the curl-type evolution, D and B satisfy involutionary constraints that are a natural part of Maxwell's equations. Maxwell's equations can be written as stress by the second rank tensor σ . The equations for a fluid with thermal conduction can be written as the heat conduction will behave asymptotically like the classical Fourier law for parabolic heat conduction. When the relaxation time is very large, the source term becomes irrelevant and the heat conduction will be described by purely hyperbolic heat waves or phonons, propagating with a characteristic speed that is called the second sound. The beauty of the above equations stems from the fact that they constitute a first order hyperbolic system with a source term that may indeed become stiff in certain limits. Therefore, all of the well-developed technologies that have been developed for solving hyperbolic PDE systems with stiff source terms can indeed be brought to bear on the numerical solution of the above PDE system. Furthermore, the solution method does not require the treatment of a parabolic sub-system, which can be computationally expensive. As already stated before, a formal asymptotic analysis of eqn. (1.3) shows that the above equations retrieve the Navier Stokes equations with the traditional Fourier law of heat conduction in the stiff limit when the relaxation time τ tends to zero. To complete our description of the above system, we also mention the constitutive relation for the thermal stress tensor Here h c denotes the hyperbolic speed of heat waves, i.e. the second sound. be specified by the accuracy of the numerical method. As a result, even for regions of the flow that should have no thermal conduction, there will indeed continue to be some small amount of thermal conduction. This affects the fidelity of the method and its results.
The fourth equation in eqn. (1.3) contains the involution and therefore deserves further attention. Because the evolution of a curl-free vector field J is only governed by the gradient term ( ) T ∇ ⋅ + J v we must pick a mimetic discretization that ensures this curl-free evolution. A good conceptual model for a curl is the altitude in a mountainous region. It does not matter which closed curve one takes in that mountainous region, as long as the curve is closed, the total change in elevation will be zero. This is the model that we keep at the back of our mind when studying this problem. The closed curve could be the edges of a rectangular or triangular mesh in 2D. For a 3D mesh, we have closed curves in all the faces of the element, whether the element is a cube or a tetrahedron. Along each of those faces, the circulation of the vector field J must be zero in all the situations where the vector field is required to evolve in a curl-free fashion. This is only guaranteed if the components of J are collocated at the edges of the mesh and ( ) is collocated at the corners of the mesh. But realize that we are solving a hyperbolic system, as a result ( ) have to be obtained consistently with multidimensional upwinding at the corners of the mesh. There already exists a 3D Riemann solver that does this (Balsara [16]). We, therefore, see our first point of analogy between involutionary PDEs that support divergence-free evolution and involutionary PDEs that support curl-free evolution. The former require 2D Riemann solvers that are invoked at the edges of the mesh. The latter require 3D Riemann solvers that are invoked at the corners of the mesh.
Let us now press on with our study of the last equation in eqn. (1.3). Let us focus on the term ( ) × ∇ × v J . When the curl is zero, it is irrelevant. However, when the curl is non-zero, it does affect the time rate of change of the component of J that is aligned with the edges of the mesh. How can we get the measure of the curl of a vector field? In three-dimensions, we can only do that by reconstructing the vector field in a three-dimensional fashion. (Likewise, of course, in two dimensions!) In other words, we need to start with the components of J in the edges that surround each volume element of the mesh and obtain from it a consistent value of J within the volume element. This should be done in a way that reflects, in some appropriate fashion, the curl that is already present in the faces of that mesh element. This is the problem of reconstructing a vector field consistent with its constraints. We met this problem already when we considered divergence-preserving reconstruction of vector fields in eqns. (1.1) and (1.2). For CED and MHD, very special forms of vector field reconstruction had to be invented that were consistent with the divergence constraint. We, therefore, see that we will have to pay special attention in this paper to curl-free and curl-preserving reconstruction of vector fields. We, therefore, see our second point of analogy between involutionary PDEs that support divergence-free evolution and involutionary PDEs that support curl-free evolution. The former requires a divergence-free reconstruction of normal components of the vector field that are collocated at the faces of the mesh elements. The latter require curl-free and curl-preserving reconstruction of tangential components of the vector field that are collocated at the edges of the mesh elements.
We have seen that eqn. (1.2) arises from taking the divergence of the first equation in eqn. (1.1). It is possible to show a similar concordance in curl constraint-preserving reconstruction. Let us, therefore, take the curl of the third equation in eqn. (1.3). Let us also make the definition ≡ ∇ × R J , where the vector field R is referred to as the Burger's vector field. It is easy to see from its very definition that 0 ∇ ⋅ = R ; i.e. the Burger's vector is divergence-free. Now let us take the curl of the fourth equation in eqn. (1.3). We get (see Peshkov et al. [50]) ( ) We see immediately that eqn. (1.4) guarantees divergence-free time-evolution for the Burger's vector field. Those who are familiar with the induction equation for MHD will also see the great parallels between eqn (1.4) evaluated in the limit where τ → ∞ and the MHD induction equation.
Notice from eqn. (1.2) that we have an equation for the time-evolution of the charge density. This means that in a DG scheme, all the moments of the charge density can be thought of as provided to us. For a finite volume (FV) scheme, of course, we have to reconstruct the higher moments of the charge density using charge densities in neighboring zones, and that is easy to do given that it is a scalar. Note that the zeroth moment of the charge density does not need to be reconstructed however because it is given by a discrete application of Gauss' law, (1.4)! Therefore, we see that when designing a DG scheme for curl constraint-preserving PDEs we will get an additional equation for the evolution of the curl of the vector field J that is of interest. For a finite volume (FV) scheme, of course, we have to reconstruct the higher moments of the Burger's vector field R . Furthermore, these higher moments must guarantee that 0 ∇ ⋅ = R . To make the analogy between divergence-preserving and curl-preserving schemes exact, the zeroth moment of the Burger's vector does not need to be reconstructed because it is always given to us by a discrete application of the definition, ≡ ∇ × R J, in the zone of interest.

I.3) Introduction : Plan of the Paper
The rest of the paper follows the ensuing plan. In Section II we show how the curlpreserving reconstruction can be carried out at all locations of a two-dimensional Cartesian mesh; this will include second to fourth order reconstructions. Section III extends these ideas to threedimensional Cartesian meshes. Those two sections also have demonstrations that the twodimensional and three-dimensional schemes are multidimensionally upwinded and, therefore, stable. Section IV presents a very efficient strategy for curl-preserving prolongation in AMR. Section V shows results of a von Neumann stability analysis of curl constraint-preserving DG-like schemes. Section VI shows some results from a model problem where the lack of curl-preserving reconstruction is shown to have obvious deleterious effects. Section VII shows results from curlpreserving prolongation, as it pertains to AMR. Section VIII presents some conclusions.

II) Curl-Preserving Reconstruction on a Two-dimensional Cartesian Mesh
It is easiest to get introduced to this subject in two dimensions, especially on a structured mesh. We consider this problem in five easy Sub-sections. In Sub-section II.1 we present a first order accurate reconstruction of a curl-preserving vector field. In Sub-section II.2 we present a second order accurate reconstruction of a curl-preserving vector field. In Sub-sections II.3 and II.4 we present third and fourth order extensions. Sub-section II.5 shows that the curl-free reconstruction, when combined with a two-dimensional Riemann solver, produces a properly upwinded numerical scheme.

II.1) Curl-Preserving Reconstruction on a Two-dimensional Cartesian Mesh at First Order
Let us consider what is entailed in a first order reconstruction. In keeping with the spirit of a first order finite volume scheme for fluid flow, it means that each edge of a rectangular/square zone has a component of the vector field along the direction of the edge. In simplest form, and for a unit square zone with extent ( ) [ ] [ ] , 1 2,1 2 1 2,1 2 x y ∈ − × − , this is shown in Fig. 1. Any rectangular zone can be mapped to such a square zone, so our results are perfectly general. Fig. 1 shows the collocation of vector components along the edges of a two-dimensional control volume. As evaluated over the edges of the square element, the discrete circulation is fully specified. (The mean value and its linear variation are shown along each edge in Fig. 1, in anticipation of a second order accurate reconstruction scheme. However, in this Sub-section we ignore the linear variation.) The reconstruction problem for a curl-free reconstruction consists of obtaining a polynomial-based vector field that is globally curl-free within this two-dimensional control volume. The reconstruction problem for a curl-preserving reconstruction consists of obtaining a polynomialbased vector field that matches the specified mean circulation in the zone.
From Fig. 1 we see that the bottom and top x-edges have x-components of the vector field that are given by 1 x V and 2 x V respectively. Likewise, the left and right y-edges have y-components of the vector field that are given by 1 y V and 2 y V respectively. A polynomial that holds over the entire unit square and matches the specified values at the edges is given by By taking the curl of the above vector field, we get ( ) We see that the curl, evaluated as a differential expression, gives back the discrete circulation of the vector field over the unit square shown in Fig. 1. If the discrete circulation is curl-free then it will evaluate to zero and our vector field in eqn. (2.1) will also be curl-free -i.e. the curl evaluated at each local point in the unit square is exactly zero. If the discrete circulation is not curl-free then the differential form in eqn. (2.2) matches the exact value of the discrete circulation at all locations of the unit square, which is reasonable. Observe too that ( ) It is also worthwhile to observe that if any three of the four components given by 1 x V , 2 x V , 2 y V and 1 y V are specified, and if we are told that the vector field is curl-free, then the fourth component is automatically satisfied. This is a small observation for now, but it will be expanded on in subsequent sections to guide us towards the problem of prolonging vector fields in curl-free (or curl-preserving) fashion on AMR meshes.

II.2) Curl-Preserving Reconstruction on a Two-dimensional Cartesian Mesh at Second Order
Now let us consider second order extensions. In the spirit of van Leer [59] and Kolgan [46], this is tantamount to endowing each of the edges with a piecewise linear variation. From Fig.  1 we see that the bottom and top x-edges have x-components of the vector field that are now endowed with undivided differences in the x-direction given by ( ) Similarly, from Fig. 1 we see that the left and right y-edges have y-components of the vector field that are now endowed with undivided differences in the y-direction given by ( ) respectively. (For a second order DG scheme, these undivided differences will indeed become evolutionary modes.) Let us say that we follow exactly the same game-plan as in eqn. (2.1) and write ( ) To see the problem with eqn. (2.3), let us take its curl. We get We see that even if the original vector field had a discrete circulation that was zero over the square shown in Fig. 1, the resulting curl evaluated at all points within the square will not be zero. This is because in general ( ) ( ) Having seen that a naïve attack on the problem yields nothing of value, let us renew our effort. We take inspiration from the divergence-free reconstruction of two-dimensional vector fields that was discussed in Sub-section III.1 of Balsara [5] and realize that when we are dealing with a constrained vector field, the components couple. In other words, the vector field is an entire entity and we cannot take the individual components as disjoint entities. Therefore, the xcomponent of the vector field will couple to the y-component of the vector field so as to preserve the constraints. Since we have already realized that curl-free and curl-preserving reconstruction are just two sides of the same coin, we focus on the former problem first. Let us write our vector field as Notice that all the terms that are needed for obtaining second order accuracy are already present in eqn. (2.5 We see that the first square bracket in the above equation still expresses the discrete circulation, which is exactly zero for a curl-free vector field. The second and third square brackets in the above equation can be made zero by setting ( ) ( ) ( ) ( ) Notice that from a finite difference point of view, the coefficients yy a and xx b are just higher order derivatives of the undivided differences, and as a result, the second order accuracy of eqn. (2.5) is not affected by the inclusion of these additional terms.
The analogies with divergence-free reconstruction in Balsara [5] are also worth drawing. In both cases, the first order term is an expression of the discrete constraint applied to the boundaries of the element. The inclusion of higher order terms requires additional coefficients to ensure that the differential form of the constraint is exactly satisfied at all locations within the element. Now that we have thoroughly discussed all the nuances of a curl-free reconstruction of a vector field, we are in a position to discuss how the idea goes over to a curl-preserving reconstruction. Recall that the fourth equation in eqn. (1.3) indeed has an evolutionary equation for the curl, see eqn. (1.4). Let us draw insights by studying the analogous situation in CED where FV and DG schemes for CED are well-developed. Consider the analogous eqn. (1.2) for CED where we also have an evolutionary equation for the charge density. The modes of the charge density, at least in a DG scheme for CED, would then be directly available to us; the same is true for a curl-preserving reconstruction. However, for a FV scheme for CED involving charge densities, Balsara et al. [22], [23] have shown that the charge density can be reconstructed to a high order approximation. We can benefit from the same insight here.
Notice that in a curl-preserving reconstruction, Fig. 1 shows us that the discrete circulation in the square zone is given by For a second order DG scheme, eqn. (1.4) would provide the time-evolving higher moments of the curl. For a FV scheme, we can reconstruct such a quantity for all zones of the two-dimensional mesh. Using neighboring elements, we can obtain a TVD-based or WENO-based piecewise linear, finite volume reconstruction of the circulation. Such a reconstruction would match the discrete circulation in the target zone. As a result, for the zone shown in Fig. 1, we can write the piecewise linear circulation as ( ) ( ) ( ) Matching eqns. (2.6) and (2.8) we get Comparing eqns. (2.7) and (2.9) we now have a ready correspondence between curl-free and curlpreserving reconstruction. Specifying one is tantamount to specifying the other.
It is also worth pointing out that eqns. (1.3) and (1.4) show us that the curl of the vector field explicitly participates in the time-update. Therefore, it is useful to provide explicit expressions not just for the vector field but also for its curl. Such expressions prove to be quite valuable for making numerical implementations.

II.3) Curl-Preserving Reconstruction on a Two-dimensional Cartesian Mesh at Third Order
We extend the results from the previous Sub-section to the third order case here. In addition to being useful for scheme design, this is useful for analytic work on DG schemes and their von Neumann stability. Notice first off that ( ) , x V x y in eqn. (2.5) has constant, x, y, xy and y 2 terms. Therefore, to become a truly third order reconstruction, it minimally needs an x 2 -dependent term, which will indeed be added along the x-edges. Similarly, ( ) , y V x y in eqn. (2.5) has constant, x, y, xy and x 2 terms. Therefore, to become a truly third order reconstruction, it minimally needs a y 2 -dependent term, which will indeed be added along the y-edges. Such a way of thinking shows us how each reconstruction of the curl constraint-preserving vector field at a certain order shows us the way to the reconstruction at the next higher order. At least on a structured mesh, where the polynomial terms can proliferate, this is the systematic strategy that one should pursue.
It is important to be emphatic about a point of detail that we develop in this paragraph. One may think that it is unreasonable to claim that the y 2 mode is present in x V x y in eqn. (2.5) because that mode comes purely from the constraint-satisfaction at second order. Similarly, one may think that it is unreasonable to claim that the x 2 mode is present in ( ) , y V x y in eqn. (2.5) because that mode also comes from constraint-satisfaction at second order. However, indeed those modes are truly present because this is the very idea behind a constrained vector field. The constraint basically tells us that if a mode in x y is needed in order to satisfy the curl-free (or curl-preserving) constraint, then it is indeed truly satisfied. It does not matter that it is satisfied by variation in the other vector component, because the curl-free (or curl-preserving) vector field is just a single entity. None of the components of the curl-free vector field are entire in themselves, they only exist as parts of a whole! An entirely analogous observation has been found to be true over and over again in divergence constraint-preserving reconstruction for MHD and CED (Balsara [5], [6], [8], Balsara et al. [22], [23]).
Let us now extend the curl-free reconstruction to third order. At the bottom and top x-edges of the square shown in Fig. 1 we now add piecewise quadratic modes that we denote by ( ) respectively. Similarly, at the left and right y-edges of the square shown in Fig. 1 we now add piecewise quadratic modes that we denote by ( ) The xyy a and xxy b are not mandatory for order property preservation, but we shall show shortly that they are needed in the construction of a DG scheme for curl constraint-preserving vector fields.
We can now write out the curl of the above vector field as ( ) As with eqn. (2.8), we can now reconstruct the discrete circulation up to and including quadratic variation over each zone, and write the result as This gives us the third order curl-free or curl-preserving reconstruction on a two-dimensional Cartesian mesh. To get a curl-free reconstruction, just set all the coefficients in eqn. (2.12) to zero.
We can now also notice that a third order accurate DG scheme which evolves all the modes of the Also notice that from a finite difference point of view, the coefficients yy a , yyy a , xx b and xxx b in eqn. (2.13) are just higher order derivatives of the undivided differences, and as a result, the third order accuracy of eqn. (2.10) is not affected by the inclusion of these additional terms.

II.4) Curl-Preserving Reconstruction on a Two-dimensional Cartesian Mesh at Fourth Order
Let us now make a fourth order extension. We use our idea of systematically thinking about the terms that are present in the third order reconstruction and using them to inform our choices at fourth order. Notice, first off, that ( ) , x V x y in eqn. (2.10) has constant, x, y, x 2 , xy, y 2 , y 3 and x 2 y terms. To that, along each x-edge, we will indeed add an x 3 -dependent term. However, to have full fourth order accurate reconstruction, we will still need an xy 2 term, which must indeed be added with a zone-centered collocation! In other words, by enriching the moments along each x-edge we simply cannot obtain a term with xy 2 variation, so we have to include it at a location where all the moments have validity, namely at the zone center. Furthermore, notice that ( ) , y V x y in eqn. (2.10) has constant, x, y, x 2 , xy, y 2 , x 3 and xy 2 terms. To that, along each y-edge, we will indeed add a y 3dependent term. However, to have full fourth order accurate reconstruction, we will still need an x 2 y term, which must indeed be added with a volume-centered collocation! As before, by enriching the moments along each y-edge we simply cannot obtain a term with x 2 y variation, so we have to include it at a location where all the moments have validity, namely at the zone center. We now see the value of our systematic, order-by-order approach because it has highlighted for us the fourth order terms that are supplied by enriching the basis space along the edges and the additional modes that have to be supplied volumetrically. (A similar subdivision occurs in divergence-free reconstruction for CED and MHD where we already know that at fourth order and beyond, many of the modes are face-centered, but some are volume-centered. There too, we realized that we could not enrich the space of spatial modes to obtain all the terms that are needed in a fourth order accurate Taylor series expansion. As a result, some of the modes had to be volume-centered; see Balsara and Käppeli [24] and Hazra et al. [44].) Let us now extend the curl-free reconstruction to fourth order. At the bottom and top xedges of the square shown in Fig. 1 we now add piecewise cubic modes that we denote by ( ) respectively. Similarly, at the left and right y-edges of the square shown in Fig. 1 x V x y will also trigger additional modes of the form x 3 y.
Likewise, the inclusion of a cubic y 3 -dependent term along each of the y-edges in ( ) , y V x y will also trigger additional modes of the form y 3 x. To compensate for the effect of these terms on the curl, some higher order polynomial terms have to be added. The analysis from the previous paragraph allows us to write the fourth order curl constraint-preserving vector field as Note that the modes xyy a and xxy b correspond to zone-centered modes that carry the 2 xy and 2 x y variation. However, those variations are calculated not to disturb the boundary variation. The terms associated with xyyyy a and xxxxy b are intended to ensure that the curl can be matched at all orders that are relevant to us. The curl now becomes As with eqn. (2.12), we can now reconstruct the discrete circulation up to and including quadratic variation over each zone and including the minimal number of cubic modes that arise in eqn. (2.15), and write the result as Equating like terms in eqns. (2.15) and (2.16) establishes most of the terms unequivocally, as follows:- So far, we have left xyy a , xxy b , xyyyy a and xxxxy b terms undetermined. They will only be fixed after we make the following consideration. The vector field in eqn. (2.14) includes modes that reside on the edges of the mesh and modes that are zone-centered. It is, therefore, interesting to ask how both kinds of modes (edge-centered and zone-centered) can be accommodated seamlessly in a DG scheme? Indeed, we get a new insight by addressing this question. Realize that the vector field in eqn. (2.14) can also be decomposed in terms of orthogonal Legendre polynomials as follows:- Now notice that owing to the mass matrix being diagonal for the expansion in eqn. (2.18), when using a zone-centered DG scheme, we are guaranteed that we need only one volumetric update equation for the time-evolution of the coefficient of ( ) x V x y and one volumetric update equation for the time-evolution of the coefficient of ( ) x z x x y y y y y available. Either way, whether we have a FV or DG scheme, the availability of these modes enables us to make two possible optimal choices for fixing the xyy a , xxy b , xyyyy a and xxxxy b terms. We discuss that in the next paragraph.
The first option will match all the curl terms exactly. It consists of setting:- However, it entails the compromise of equating xyy a and xxy b in the limit where ( ) 0 z xy R ∆ → . In other words, the coefficient of ( ) The second option will not force the coefficients to coincide. But it will match all the curl terms exactly, up to the desired order of accuracy while minimizing (in a least squares sense) the curl terms that are not mandated by order of accuracy considerations. This choice consists of setting:- For the test problems that we have considered in the results sections, both these choices have worked well and retained the order of accuracy. We will show results for the latter choice in eqn. (2.21). This completes our description of curl-preserving reconstruction at fourth order in two dimensions.

II.5) Combining Curl-Free Reconstruction and the Two-Dimensional Riemann Solver to Obtain a Multidimensionally Upwinded, Globally Curl-Free Scheme
When studying one-dimensional advection, it is indeed a very instructive to realize that one-dimensional upwinding from a one-dimensional Riemann solver yields a stable parabolized scheme. The transition to a higher order scheme for advection is then easy to justify with the inclusion of TVD or WENO limiters, as was shown by van Leer [59], [60], Jiang and Shu [45], Balsara and Shu [4], Balsara, Garain and Shu [19]. For third order accurate reconstruction, the PPM schemes of Colella and Woodward [33], Colella and Sekora [34] or McCorquodale and Colella [48] would work as well. As the order of accuracy of the reconstruction is increased, it is easy to see that the dissipative terms will become smaller. It is very desirable to show that an analogous plan exists for PDEs with an involution constraint on a vector field. Such a demonstration has to be multidimensional because the curl operator is only meaningful in two or more dimensions.
In their study of the induction equation for globally divergence-free MHD, Balsara and Käppeli [20] were able to show that multidimensional divergence-free reconstruction of the magnetic field at first order, coupled with a two-dimensional Riemann solver, also results in a stable scheme with the correct, multidimensionally parabolized, dissipation; see Section 4 of the above-mentioned paper. Because of the presence of these parabolic terms, the multidimensional Riemann solver always plays a stabilizing role in the induction equation. This showed us that the transition to higher order, by applying TVD, PPM or WENO-type limiting, would indeed succeed for the induction equation. Notice that Balsara and Käppeli [20] have made a non-trivial demonstration, because the components of a divergence-free vector field are collocated at the faces of the mesh, which is quite different from the zone-centered collocation that is used in DG schemes for conservation laws.
The previous paragraph has shown us that divergence-free evolution of vector fields has been developed to a point where it has a solid footing. It is very desirable to show that vector fields that evolve in a curl-free fashion have a similar assurance. This will again be a non-trivial demonstration because the components of a curl-free vector field are indeed collocated at the edges of the mesh. This is indeed different from the zone-centered collocation that is used for conservation laws as well as the face-centered collocation that is used for mimetic schemes that support globally divergence-free evolution. (Curl constraint-preserving evolution is then just a matter of adding source terms to the curl-free evolution equations; so we do not consider that in this Sub-section.) Let us consider the simplest two-dimensional equations that give us curl-free evolution. They can be written as Here we will take the velocity components ( ) v , v x y to be constant. Let us view the update equations shown above in two space and one time direction. We also consider a uniform Cartesian mesh in the two spatial dimensions with zones of size x ∆ and y ∆ in the x-and y-directions. Let the timestep be of size t ∆ . The first equation in eqn. (2.22) can be integrated along the x-and tdirections of the three-dimensional space-time mesh (2 space + 1 time dimensions). Since x J is collocated at the x-edges of the mesh, the x-directional integration should, of course, coincide with the x-edges of the Cartesian mesh. As long as the curl is exactly zero, the term in that equation will be exactly zero; and the ( ) applied in the x-direction. Similarly, the second equation in eqn. (2.22) can be integrated along the y-and t-directions of a three-dimensional space-time mesh. Since y J is collocated at the y-edges of the mesh, the y-directional integration should, of course, coincide with the y-edges of the two-dimensional Cartesian mesh. As long as the curl is exactly zero, the ( ) , are used at the corners of the mesh in order to update the x-edge centered x J and the y-edge centered y J then the update will be globally curl-free. (This is very similar to globally divergence-free update of the induction equation in MHD which requires the same, properly upwinded, edge-centered, electric fields to be used for the update of the face-centered components of the magnetic induction vector.) Let us now consider constant velocity components with v 0 x > and v 0 y > , just to keep our initial discussion simple. Eqn.
We see, therefore, that "properly upwinded" in this context means that each zone in Fig. 2 will contribute its x J and y J values to its top right corner in Fig. 2. We require that the first order accurate curl-free reconstruction from eqn. (2.1) should be used in the zones ( ) J − indicate spatial collocation points on the mesh; the superscript of "n" denotes the n th timestep. The components of the vector field J are also indicated by "x" or "y" in the superscripts. We can write the discretized, first order in space and time, update equations as Recall that we assume a uniform Cartesian mesh with zone sizes x ∆ and y ∆ in the x-and ydirections and a timestep of size t ∆ . By substituting the upwinded potentials from Fig. 2 We now write the above equations in a format that allows us to see the velocity-dependence more clearly as follows The above equations still do not look like discretized versions of eqn. (2.23). As with the divergence-free evolution of vector fields, the concordance will only be established if the discrete circulations in the zones ( ) are used in the first and second equations of eqn.
(2.26). If we utilize the fact that the discrete circulations in those two zones are indeed zero, we can make the transcriptions It is now easy to see that the first equation in eqn. The two innovations work together to yield the globally curl-preserving scheme. We therefore see with the help of Fig. 2 that the curl-free reconstruction, along with multidimensional upwinding applied to the vertices of the mesh, gives us a globally curl-free update strategy for our model equations, i.e. eqn. (2.22). Realize too that the multidimensional Riemann solver is indeed designed to automate that multidimensional upwinding in the general case of a system of PDEs. Therefore, we realize that the curl-free reconstruction (and the edge-centered collocation of vector components that it entails), along with the application of a multidimensional Riemann solver, indeed gives us a stable, globally curl-free, mimetic update strategy for our curl-free model equations.
The discussion in the previous paragraph was restricted to first order accuracy. To keep the discussion extremely accessible, we also drew on our notional understanding of multidimensional upwinding, and we kept all the velocities positive. In general, we will want to use a full multidimensional Riemann solver which can accommodate to waves propagating in a general hyperbolic system in any direction. Eqns. (12), (13) and (14) of Balsara [13] show how such a multidimensional Riemann solver can be designed for systems in the general case. The multidimensional Riemann solver from Balsara [13] works on Cartesian meshes. It was extended to unstructured meshes in Balsara and Dumbser [15]. For our model equations in eqn. (2.22), we do not have the support of a general underlying PDE system, however, we can design a multidimensional Riemann solver (using eqns. (12), (13) and (14) of Balsara [13]) in the locally Lax-Friedrichs (LLF) approximation. This enables us to write our potentials at the vertices of Fig.  2 as a centered part as well as a multidimensionally dissipative part. We will then write discrete evolution equations and show that we get a centered update along with a parabolic contribution at first order. This makes it easier for us to understand that as the order of accuracy of our curl-free reconstruction is improved, and as the accuracy of our time-update is improved, the dissipative parts will become progressively smaller. In other words, the same exercise that was presented in Section 4 of Balsara and Käppeli [20] for mimetic, globally divergence-free schemes is now replicated in the ensuing paragraphs for mimetic, globally curl-free schemes.
We still consider the case where the velocity components ( ) v , v x y are constant but now they can have any sign. We write the potential at the vertex ( ) Analogous expressions can be written for 1 The above equations clearly show us the centered terms and the dissipation terms in the first order scheme. We see that the dissipation terms have the correct scaling to perfectly stabilize the scheme. The use of higher order curl-free reconstruction and higher order time stepping, in conjunction with the multidimensional Riemann solver, will reduce the dissipation. We now have our assurance that higher order, mimetic, globally curl-free schemes will be stable.

III) Curl-Preserving Reconstruction on a Three-dimensional Cartesian Mesh
The prior exercise at reconstructing curl-free or curl-preserving vector fields in two dimensions has left us with many very valuable insights. We list them below:-

1)
We see that one should start at the lowest order and systematically build up to higher orders. This is because the polynomial terms at each order give us insight into which modes are needed at the next higher order.

2)
We also obtained the insight that some additional polynomial contributions will be needed to ensure curl-free vector fields at all locations within a zone. Because of the nature of the curl operator, and our need to only use the lowest order additional polynomials so as to retain a modicum of stability, we found that one component usually takes on contributions that help cancel extra terms that arise in another component.

3)
At fourth and higher orders the structure of the mathematics is such that it requires some modes to be evolved in zone-centered fashion when designing a DG scheme. In such situations, the presence of a diagonal mass matrix is very helpful.
Armed with these insights, we now extend our studies to three-dimensional Cartesian meshes.
Sub-section III.1 deals with the first order reconstruction on 3D Cartesian meshes. Subsection III.2 extends this to second order. Sub-sections III.3 and III.4 present the third and fourth order cases. Sub-section III.5 shows that the curl-free reconstruction, when combined with a threedimensional Riemann solver, produces a properly upwinded numerical scheme. Within each of the two x-faces, the two y-faces and the two z-faces, the discrete circulation (evaluated over those faces) is exactly zero if we are dealing with a curl-free PDE. If the PDE is only curl-preserving, the discrete circulation is easily evaluated at each face, as we soon show. The mean value of the vector field's parallel component and its linear variation are shown along each edge, in anticipation of a second order accurate reconstruction scheme. The reconstruction problem consists of obtaining a polynomial that is globally curl-free/curl-preserving within this control volume. The discrete circulation within each of the six faces of the zone shown in Fig. 3 gives us the six conditions for the mean values. For simplicity, let us consider curl-free evolution in the next three equations. At the top and bottom x-faces we have

III.1) Curl-Preserving Reconstruction on a Three-dimensional Cartesian Mesh at First Order
At the top and bottom y-faces we have At the top and bottom z-faces we have In general, the second and higher moments in Fig. 3 are only obtained with some level of approximation. So we cannot guarantee that an analogous set of equations hold for the slopes.
It is very important to make one further observation about a curl-free vector field. The six discrete curl conditions in eqns. (3.1) and can be written as 3 y V appears only in the second discrete circulation condition in eqn. (3.1) and can be written as x V only appears in the first discrete circulation condition in eqn. (3.2) and can be written as x V only appears in the second discrete circulation condition in eqn. (3.2) and can be written as 1 Adding the above four equations immediately retrieves the first discrete circulation condition in eqn. (3.3). We have started with the second discrete circulation condition in eqn.  .2) can be assumed true. Therefore, on a mesh with six faces, only five of them are truly, mutually independent. We will see that this plays a very important role in AMR, as will be shown in Section IV.
Consider the following reconstruction on the unit cube spanning [ ] We evaluate its curl in the next paragraph.
in the above equation to see that it retrieves the discrete circulation equations in eqn. (3.1) at the top and bottom x-faces respectively. If the vector field is circulation-free then this in the above equation to see that it retrieves the discrete circulation equations in eqn.
in the above equation to see that it retrieves the discrete circulation equations in eqn. (3.3) at the top and bottom z-faces respectively. If the vector field is circulation-free then this guarantees that ( ) z ∇ × V is strictly zero in the above equation.
(3.7), (3.8) and (3.9) we also see that In other words, the discrete divergence of the discrete curl is also exactly zero. This also tells us that a good zone-centered approximation of ( ) x ∇ × V is given by the first term in eqn. (3.7). Likewise, a good zone-centered approximation of ( ) y ∇ × V is given by the first term in eqn. (3.8). Similarly, a good zone-centered approximation of ( ) z ∇ × V is given by the first term in eqn. (3.9). This tells us that if we want to build higher order approximations of the curl, we could start with the identified terms in each zone and then use higher order WENO reconstruction to obtain the higher moments. However, the slight nuance is as follows. If the higher modes of a reconstructed curl operator are obtained via a finite volume approximation, we should also make sure that they are divergence-free.
Now notice that eqns. (3.4), (3.5) and (3.6) are only first order accurate. This is because eqn. (3.4) lacks any linear variation in the x-direction; eqn. (3.5) lacks any linear variation in the y-direction and eqn. (3.6) lacks any linear variation in the z-direction. However, for the first order accurate case, the equations are exactly curl-free or curl-preserving. Besides, the first order curlfree reconstruction reflects the six discrete circulations evaluated over the six faces of the control volume using the edges of the same control volume. (This mirrors the known fact that at first order, the discrete divergence-preserving reconstruction reflects the one discrete divergence evaluated over the control volume using the faces of the same control volume.) Notice too that while there is only one divergence condition evaluated over a 3D control volume in a divergence-preserving scheme, there are six curl conditions evaluated over a 3D control volume in a curl-preserving scheme. This makes curl-free reconstruction more complicated than divergence-free reconstruction, especially in three dimensions.
Eqns. (3.7), (3.8) and (3.9) give us yet another insight if we compare them to eqn. (2.2). Notice that eqn. (2.2) is a single equation for the discrete circulation. However, because of the linear variation in the x-direction, eqn. (3.7) is an expression of the discrete circulation at the top and bottom x-faces. Similarly, because of the linear variation in the y-direction, eqn. (3.8) is an expression of the discrete circulation at the top and bottom y-faces. Likewise, because of the linear variation in the z-direction, eqn. (3.9) is an expression of the discrete circulation at the top and bottom z-faces. We, therefore, understand that even when we make a higher order reconstruction of the curl vector in the zone of interest, the two modes that are present in each of eqns. (3.7), (3.8) and (3.9) must be retained.

III.2) Curl-Preserving Reconstruction on a Three-dimensional Cartesian Mesh at Second Order
Notice that eqn. (3.4) already has a constant part and y, z and yz variation. To attain full second order accuracy, we need to add a linear x-directional variation to ( ) , , x V x y z . This would be added to the x-edges of the zone shown in Fig. 3. The inclusion of such an x-variation will also trigger xy, xz and xyz terms in ( ) , , x V x y z . Similarly, notice that eqn. (3.5) already has a constant part and x, z and xz variation. To attain full second order accuracy, we need to add a linear ydirectional variation to ( ) , , y V x y z . This would be added to the y-edges of the zone shown in Fig.   3. The inclusion of such a y-variation will also trigger xy, yz and xyz terms in Likewise, notice that eqn. (3.6) already has a constant part and x, y and xy variation. To attain full second order accuracy, we need to add a linear z-directional variation to ( ) , , z V x y z . This would be added to the z-edges of the zone shown in Fig. 3. The inclusion of such a z-variation will also trigger xz, yz and xyz terms in The inclusion of all these terms also causes the curl operator to acquire additional moments and to ensure curl-free reconstruction (or to make sure that the curl-preserving reconstruction has the appropriate moments) we need to add some complementing terms.
Notice too that the polynomials are chosen with such factors that they do not affect the vector components in the edges. In other words, if the ensuing polynomial for    Now let us study the curl-preserving constraints that are put on the above vector field. Appendix A catalogues the curl of this second order accurate vector field along with some helpful discussion on how the terms are to be matched.
Because the second order accurate vector field only has to be specified up to linear terms, its discrete curl only needs to be specified up to constant terms when we are considering a second order finite volume scheme. But realize from an examination of eqn. (1.4) and how it arises from the last equation in eqn. (1.3) that a second order DG scheme will also have an evolutionary equation for the curl of the original vector field; and the latter equation also needs to be evolved with second order accuracy. (In two-dimensions, this consideration was immaterial; but in three dimensions it is a meaningful consideration.) However, for a DG scheme, the curl vector itself might have been evolved up to piecewise linear terms as follows Note that the first two terms in the above three equations are mandatory and cannot be changed. This is because they are needed to retrieve the discrete circulation aspect of the curl vector and also to ensure that the curl vector is divergence-free. (This is analogous to the charge density in any FV/DG scheme for CED. The mean value of the charge density in any zone must indeed be compatible with the discrete divergence of the electric displacement evaluated from the faces of the mesh.) We now equate eqns.
We see immediately that the serendipitous cancellation that we sought is indeed achieved. We also get It is easy to see the finite difference-like structure in the above nine equations. We see that they truly represent higher derivatives that can be derived from the arrangement of gradients at the edges of the control volume in Fig. 3. This also tells us that the effect of these higher derivative terms in eqns. (3.10), (3.11) and (3.12) is to slightly modify the higher moments in those equations so as to restore curl-free or curl-preserving behavior. It is also important to note that, in spite of this modification, the reconstructed vector field will indeed exactly match the values of the gradients at the edges. Note too that, the modifications will be slight owing to the fact that the modifying coefficients represent higher derivatives. Therefore, the reconstructed vector field in eqns. (3.10), (3.11) and (3.12) is suitable for the construction of curl constraint-preserving schemes. This completes our study of curl constraint-preserving reconstruction of vector fields at second order of accuracy.
For second order accurate AMR with a refinement ratio of two, the reconstruction given in this Sub-section can be directly applied to any unrefined zone that needs to have its data prolonged to a finer mesh in a curl-preserving fashion. However, it is also possible to design a curl constraintpreserving AMR strategy based on any sufficiently accurate reconstruction, including one that does not satisfy the curl constraint (to begin with). The approach in Section IV works as long as the reconstruction we start with is high order accurate, in other words, the starting reconstruction does not have to be curl constraint-preserving and may be obtained from a simple finite volume WENO of the desired accuracy. Using the innovative touch-up procedure in Section IV, we can extend it to any AMR refinement ratio and up to the order of the starting reconstruction. Furthermore, we can do this in a way where the final prolongation will nevertheless be curlconstraint preserving.

III.3) Curl-Preserving Reconstruction on a Three-dimensional Cartesian Mesh at Third Order
To understand the nuances that enter the reconstruction procedure at third order, it is helpful to divide our study into two stages in the following two paragraphs. In the first stage, it is helpful to study curl-free reconstruction of vector fields and the restrictions it places on the various terms in the reconstruction. Only in the second stage will we study curl constraint-preserving reconstruction of vector fields. This two-stage sub-division is also useful because several very useful PDE systems, like general relativity, only require a curl-free reconstruction of vector fields and do not need the computationally heavier details of a curl-preserving reconstruction of vector fields.
In this first stage of our study, let us begin by considering a curl-free vector field that is represented on a three-dimensional Cartesian mesh with third order of accuracy. The component ( ) , , x V x y z in eqn. (3.10) has the following modes:-It has a constant mode and it also has modes for   2  2  2  2 , , , , , , , , , , x y z xy xz yz y z y z yz xyz . Notice that the only mode that is missing in ( ) , , x V x y z for obtaining full third order accuracy is the 2 x mode. This is good news, because it means that along the four x-edges of Fig. 3 we have to reconstruct/evolve the second moment of the x-component of the vector field along those edges. Let us label those four modes ( ) with an obvious extension of the notation. (This extension of notation also applies to the y-and z-directions. ) When an 2 x mode is added along each of the four x-edges, it also triggers the additional presence of 2 2 2 , , x y x z x yz variation in ( ) , , x V x y z . Similar considerations apply to the other directions so that we realize that along the four y-edges of Fig. 3 we have to provide four 2 y -dependent modes to In this second stage of our study, we realize that the evolutionary equations for a curlpreserving vector field will have a corresponding evolutionary equation for its curl vector. For instance, see eqn.  .14) and (3.15); and furthermore, eqn. (3.17) shows us that these additional modes have an impact on the reconstruction.) We realize that at third order we will similarly have to include additional modes in the curl-preserving reconstruction of the vector field. However, the number of such modes that we have to include is quite large. For that reason, we draw a distinction between the minimum number of modes that are needed for curl-free reconstruction and the modes that are additionally needed for curl-preserving reconstruction. That way, a user who is implementing just a curl-free reconstruction does not have to include the additional modes. To make this distinction very clear we make the following convention:-The terms which begin with the small letters "a", "b" or "c" are needed for curl-free reconstruction. The terms which begin with the capital letters "A", "B" or "C" are extra terms that are only needed for curl-preserving reconstruction.
As in the second order case, the polynomials are chosen with such factors that they do not affect the vectors in the edges. In other words, if the ensuing polynomial for     Now let us study the curl-preserving constraints that are imposed on the above vector field. Appendix A catalogues the curl of the above-mentioned third order accurate vector field.
Because the third order accurate vector field only has to be specified up to quadratic terms, its discrete curl only needs to be specified up to linear terms when we are considering a third order finite volume scheme. But realize from an examination of eqn. (1.4) and how it arises from the last equation in eqn. (1.3) that a third order DG scheme will also have an evolutionary equation for the curl of the original vector field; and the latter equation also needs to be evolved with third order accuracy. However, for a DG scheme, the curl vector itself might have been evolved up to piecewise quadratic terms as follows Notice that the ( ) 6 x xx R − ∆ in the constant part of ( ) , , x R x y z is needed in order to obtain the correct value of the discrete circulation at the top and bottom x-faces of the zone in Fig. 3. Please note that in a DG scheme, we can evolve the higher moments in the above equation. However the constant term as well as the linear in "x" term in the above equation are not evolved. They are reset using the edge values of the vector field at the end of every timestep. We also have As before, the ( ) 6 y yy R − ∆ in the constant part of ( ) , , y R x y z is needed in order to obtain the correct value of the discrete circulation at the top and bottom y-faces of the zone in Fig. 3. Please note that in a DG scheme, we can evolve the higher moments in the above equation. However the constant term as well as the linear in "y" term in the above equation are not evolved. They are reset using the edge values of the vector field at the end of every timestep. We also have As before, the ( ) 6 z zz R − ∆ in the constant part of ( ) , , z R x y z is needed in order to obtain the correct value of the discrete circulation at the top and bottom z-faces of the zone in Fig. 3. Please note that in a DG scheme, we can evolve the higher moments in the above equation. However the constant term as well as the linear in "z" term in the above equation are not evolved. They are reset using the edge values of the vector field at the end of every timestep. It is also important to realize that the vector field given in the above three equations is divergence-free. As a result, we also have the divergence-free constraints There is a further set of terms which are needed to satisfy the highest power terms in the curl constraint equations. They also depend on the coefficients in eqns.
There is a further set of terms which indeed depend on the coefficients in eqns. ( The extra terms that we introduced in eqns. ( This completes our study of curl constraint-preserving reconstruction of vector fields at third order of accuracy.
For third order accurate AMR with a refinement ratio of two, the reconstruction given in this Sub-section can be directly applied to any unrefined zone that needs to have its data prolonged to a finer mesh in a curl-preserving fashion. However, it is also possible to design a curl constraintpreserving AMR strategy based on any sufficiently accurate reconstruction, which works with any refinement ratio, as shown in Section IV.

III.4) Curl-Preserving Reconstruction on a Three-dimensional Cartesian Mesh at Fourth Order
At fourth order, and in three dimensions, the analysis becomes very tedious. In the next section we will discuss the constraint-preserving prolongation problem in AMR. There we will establish an important connection between the prolongation problem and the reconstruction problem. It will provide an alternative process of achieving constraint-preserving reconstruction. The method described there only requires that we have some procedure for starting with the edgecentered variables and their higher moments and using them to obtain a volume-averaged value for the vector field. The strategy for doing that was already shown in two-dimensions and at fourth order in eqn. (2.19), and analogous expressions for three-dimensions are given below:- 24 We see that the above expressions are easy to evaluate. We just point out that the terms like ( ) etc. in the above expressions should be obtained with at least third order accuracy (i.e. the reconstruction should be a fourth order accurate reconstruction) if the fourth order property is to be retained. Once the above three volume averages are obtained, Section IV provides a very simple strategy for the curl-preserving reconstruction that builds on just a higher order finite volume reconstruction. Since multidimensional, fourth order, finite volume reconstruction is a technology that everyone has mastered, we show in Section IV that a small tweak to this well-known technology can give us higher order curl-preserving reconstruction.

III.5) Combining Curl-Free Reconstruction and the Three-Dimensional Riemann Solver to Obtain a Multidimensionally Upwinded, Globally Curl-Free Scheme
By now it is very evident that curl-free reconstruction of vector fields requires a collocation of vector components at the edges of the mesh. This is true in two and three dimensions. In Subsection II.5 we showed that in two dimensions, the use of curl-free reconstruction, in conjunction with a two-dimensional Riemann solver that provides the two-dimensional upwinding, can indeed result in a properly upwinded, globally curl-free scheme. Since the curious reader might wonder whether there is an analogous extension to three dimensions, we provide such an extension here. To keep the discussion generally applicable to any possible time stepping strategy, we present it within the context of a semi-discrete formulation in time.
Let us focus on the model system that is the three dimensional extension of the one we studied in eqn. (2.22). We have Analogous to the demonstration in Sub-section II.5, it is easy to show that if a curl-free vector field is initially provided, and if a curl-free reconstruction is used, then the semi-discrete form of the above equations ought to evolve the vector field in curl-free fashion. Additionally, we require that the potential be collocated at the vertices of the mesh. The above system is of great practical interest because it arises naturally as part of the first order CCZ4 hyperbolic system that has to be solved for numerical general relativity; in fact, in the general relativistic context the equations do not have source terms. As in Sub-section II.5, the emphasis in this Sub-section is on curl-free evolution of the three-dimensional vector field because the inclusion of source terms on the right hand sides of the above equations can easily turn them into curl-preserving equations.
While the above equations suggest that we are simply advecting a vector field in three dimensions, please note that eqn. (3.30) enjoins us to keep the vector field globally curl-free. This requires the edge-centered collocation of the components of the vector field. To obtain stable advection, we need to show that the update methodology is also properly upwinded. We demonstrate all these facets within the context of a first order semi-discrete scheme. Fig. 4 shows the collocation of curl-free vector components along the edges of a threedimensional zone. The zone center is indexed by (i,j,k) and the edges are indexed suitably, consistent with the zone center's indexing. As in the two-dimensional case, the potentials are defined by and they are collocated at the vertices of the mesh. As long as the same potential at a vertex is used for the update of all the vector components in all the edges that meet at that vertex, the update will be globally curl-free. To keep the discussion simple, we take all the velocity components to be constant and positive. All the zones of the Cartesian mesh are also taken to be uniform with mesh sizes x ∆ , y ∆ and z ∆ in the x-, y-and z-directions. The upwinded potentials at two of the vertices of the mesh are also shown. The potentials at other vertices can be obtained by suitable shifts in the indexing. The purpose of this figure is to make it easy for us to understand how a curl-free reconstruction that is based on edge-centered vector components, in conjunction with a three dimensional Riemann solver, can give us a stable, globally curl-free scheme.
In the ensuing discussion, we shall focus on only the first of the equations in eqns. (3.30) and (3.31) because manipulations that are identical to the ones shown below can be made for the other two components of the vector field J . Since the vector field is curl-free, the terms The concordance between eqn. (3.35)

and the first equation in eqn. (3.31) is now very obvious.
We see that at first order our globally curl-free scheme is indeed properly upwinded. For the update in eqn. (3.35) we should use an effective control volume of size x y z ∆ × ∆ × ∆ that is centered on the center of the edge ( ) , 1/ 2, 1/ 2 i j k + + shown in Fig. 4. As we use higher order curl-free or curl-preserving reconstructions and higher order timestepping, we have the assurance that the numerical dissipation will be progressively reduced with increasing order but the scheme will remain stable. The demonstration might require a few steps, but the previous equation can also be written in a form that makes the centered and dissipation parts self-evident.
We now see that a higher order extension needs three essential ingredients. We need:-1) a higher order curl-free reconstruction of the sort that is presented here, 2) the three-dimensional Riemann solver from Balsara [16] and 3) a higher order timestepping strategy. This combination of innovations will produce a stable, high order, multidimensionally upwinded, globally curl-free scheme.

IV) The Role of Constraint-Preserving Reconstruction in the Prolongation Problem for AMR
Adaptive Mesh Refinement (Berger and Colella [30]) is a very powerful technique for solving PDEs for problems on a nested hierarchy of adaptive meshes where the interesting physics has to be captured on a range of length scales and/or a range of time scales. Early AMR methods were developed for fluid flow where the mass, momentum and energy density of the fluid were the primal variables. Because these are volumetric densities, it was realized that the same finite volume reconstruction methods that are used for reconstructing the fluid variables with high orders of accuracy in a higher order Godunov scheme could indeed be used for the prolongation of flow variables from a coarse mesh to a finer mesh in an AMR hierarchy. While there are other steps in an AMR scheme (such as flux correction and timestep sub-cycling), the prolongation problem tends to be the hardest one to solve. This became apparent when AMR-MHD was introduced in Balsara [5]. The magnetic field preserves a divergence-free constraint and this constraint in the vector field has to be preserved as one prolongs the magnetic field from coarser meshes to finer meshes in an AMR hierarchy. Therefore, while there are other steps in an AMR-MHD scheme (such as electric field correction and timestep sub-cycling), the high order divergence constraintpreserving prolongation problem was again the hardest problem that one had to solve. We, therefore, anticipate that the curl-preserving prolongation of vector fields will again be the hardest problem that one has to solve in AMR for such curl-constrained PDE systems.
Our goal in this section is to show that the solution to the problem of prolonging vector fields in a curl-preserving fashion can be built from the simplest of considerations. To this end, realize that higher order finite volume reconstruction is a technology that has been mastered by practically all research groups. We show here that on a coarse mesh we can start with the edgecentered primal variables that are central to curl constraint-preserving schemes and obtain from them the corresponding volume-averaged vector fields. To get the exact expressions for the volume-averaged vector fields, please see eqn. (2.19) for the two-dimensional case and eqn. (3.29) for the three-dimensional case. Both these equations are fourth order accurate as long as the onedimensional moments in the edges are constructed with sufficiently high order of accuracy. This is accomplished with a high order one-dimensional reconstruction along each edge. Starting with these volume averages within each zone, we perform a three-dimensional finite volume reconstruction at high enough order, say using a WENO or higher order PPM method. In this work we used WENO reconstruction. The resulting vector field can then be interpolated to all locations on the fine mesh. At the fine mesh edges that coincide with the coarse mesh edges, we indeed use the edge-centered one-dimensional reconstruction on the coarse mesh to provide the best possible values to the fine mesh. However, we will show that at other locations, the higher order, threedimensional finite volume WENO reconstruction is sufficient. For shared faces on the coarse mesh, one can arithmetically average from both sides of a face. For the edges of the fine mesh that do not coincide with the edges or faces of the coarse mesh, just the FV reconstruction can be used. This gives us a vector field on the fine mesh that is not curl-preserving. However, the resulting field has some very interesting properties. We study that vector field via an example below because it guides us to an extremely efficient way to obtain curl-preserving prolongation on the fine mesh by performing a small extra touch-up step to the prolonged fine mesh vector field. We set up this numerical example in the paragraphs below. The next couple of subsections will detail this touchup step in two and three dimensions.
Let us illustrate the viability of a touch-up step with a numerical example. We start our numerical example by considering the scalar field on the unit cube given by x y z x y z y z x z x y ϕ π π π π π π π π π = + − + − + (4.1) We first assign this field to the vertices of a coarse mesh that spans the unit cube. Then we difference it along each edge to find the vector field along each edge of the coarse mesh that is exactly curl-free. We now use one-dimensional WENO reconstruction along each of the three directions to build a high order reconstruction of the component of the vector field along each edge. Using eqn. (3.29) then allows us to build a zone-averaged representation of the vector field with sufficiently high order of accuracy. (We restrict ourselves to second, third and fourth orders of accuracy here.) A very standard three-dimensional finite-volume (FV) WENO reconstruction with the desired order of accuracy is now applied to the each of the three components of this vector field.
Now realize that the fine mesh has three types of edge locations. First, there are the edges of the fine mesh that coincide with the edges of the coarse mesh. The one-dimensional WENO reconstruction along each coarse edge can be used to make an exact assignment of the component of the vector field to the coincident edges of the fine mesh. Second, there are edges on the fine mesh that lie within the faces of the coarse mesh. For those fine mesh edges, we can arithmetically average the appropriate components of the 3D FV reconstruction from either of the two coarse mesh zones that share the coarse mesh face. For example, there will be x-directional and ydirectional edges from the fine mesh that lie within the xy-face of the coarse mesh. Those x-edges and y-edges will get their x-and y-directional vector fields respectively from an arithmetic averaging process applied to the three-dimensionally reconstructed vector field within the two coarse mesh zones that come together at that coarse mesh face. Third, there are edges of the fine mesh that lie entirely within a coarse mesh zone. These edges will get their vector field components from the volumetrically reconstructed vector field within the parent coarse mesh.
Notice that because the FV WENO reconstruction does not preserve the curl-free constraint, we do not expect the resulting vector field on the fine mesh to be curl-free. However, let us look at the results of the numerical experiment. In this numerical experiment, we initialize the vector field that results from eqn. (4.1) on a coarse mesh and then we prolong that vector field to the fine mesh. We use the simple procedure described in the previous paragraph for the prolongation. Because the prolonged vector field can be directly compared with the analytic answer evaluated on the fine mesh, we can find the order of accuracy of our prolongation process. We can also find the magnitude of the circulation within each face of the fine mesh. The discrete circulation evaluated on the fine mesh will indeed be non-zero over the faces of the fine mesh, but again the point of interest is to measure the order of accuracy with which the discrete circulation shrinks with increasing mesh refinement. We observe a very interesting trend. Table I shows the L1 and L ∞ errors in the x-component of the prolonged vector field as well as the maximum absolute value of the circulation on the mesh. Table I was produced by using second, third and fourth order FV WENO reconstruction with a refinement ratio of two. The WENO schemes we used were very standard schemes of finite volume type (Balsara, Garain and Shu [19], Balsara et al. [7]) and nothing special was done with respect to the curl. We see that the prolonged vector field on the fine mesh meets the appropriate order of accuracy. However, we also obtain the very important insight that the discrete circulation on the fine mesh also meets the appropriate order of accuracy! Moreover, we see that the discrete circulation, while non-zero, diminishes very quickly with mesh refinement. This gives us the insight that if a small touch-up procedure can be invented to mildly redistribute the vector components within the edges of the fine mesh by just a little then we can make the solution globally curl-free on the fine mesh.
This touch-up redistribution should be done in an optimal way, so as to refrain from destroying the order of accuracy. This redistribution should also be completely local so that it does not involve the inversion of global matrices. We are helped in that effort when we realize that the first type of edges of the fine mesh, i.e. the ones that coincide with the edges of the coarse mesh, should not have their values changed. As a result, all changes should be localized to the second type of fine mesh edges or the third type of fine mesh edges. Recall that the second type of fine mesh edges lie entirely within the faces of the coarse mesh and the third type of fine mesh edges lie entirely within the zonal volumes of the coarse mesh. Each set of 4 second type of fine mesh edges within a face is bounded by 8 type one fine mesh edges; please preview Fig. 5. Each set of 6 third type of fine mesh edges within a volume is bounded by 24 type one fine mesh edges and 24 type two fine mesh edges; please preview Fig. 7. Table I shows the result of using a very standard finite volume WENO reconstruction at second, third and fourth order to prolong a curl-free vector field from a coarse mesh to a finer mesh with a refinement ratio of two. Only the x-component of the curl-free vector field is shown because the other components show a similar trend. We show the errors in the prolongation of the vector field; furthermore, we also show the maximum circulation of the vector field on the faces of the fine mesh. Because the FV WENO reconstruction is not curlfree, we do not expect the prolonged vector field to the curl-free on the finer mesh. Sub-section IV.1 shows how this touch-up procedure can be applied to the fine mesh edges of second type. It restores curl-free/curl-preserving structure in two-dimensions. Sub-section IV.2 shows how this touch-up procedure can be applied to the fine mesh edges of third type. It restores curl-free/curl-preserving structure in three-dimensions.

IV.1) Restoring Curl-Preserving Structure in Two-Dimensions
We now describe the strategy for carrying out a curl-preserving prolongation of a vector field on a two-dimensional mesh. (It is simplest to understand the two-dimensional case before moving on to the three-dimensional case.) The same logic is also used for prolonging fine mesh edges of type two, i.e. edges that reside within the face of a coarse zone, for three-dimensional prolongation. While the method we present is instantiated for a refinement ratio of two in the narrative, towards the end of this Sub-section, we will generalize it to accommodate any refinement ratio. We will also show how this prolongation strategy yields a very general curl-preserving reconstruction strategy that can be extended to all orders.
Please focus on Fig. 5. In Fig 5, the edge numbering is loosely patterned after Fig. 1. The solid black lines show one zone of a coarse two-dimensional mesh. The dashed magenta lines demarcate the four fine mesh zones that are obtained by sub-dividing the coarse mesh zone with a refinement ratio of 2. The fine mesh zones are labeled 1 to 4. At each edge of the coarse mesh, a one-dimensional, suitably higher order WENO reconstruction is applied to the vector field component in the direction of that edge. This gives us the reconstructed component of the vector field along that edge. As a result, the one-dimensional reconstruction along each edge of the coarse mesh can be used to obtain two higher order interpolated values along that edge. These are shown as 11 These values remain the same across adjoining zones of the fine mesh; giving us reason to think that a globally curl-preserving touchup strategy can be invented for the fine mesh in two dimensions. Because of this concordance across adjoining zones, we do not have to resort to a global matrix inversion when enforcing the constraint. The constraint is only enforced locally within the four internal fine mesh edges that are formed by subdividing the coarse mesh.
Also realize that because of the higher order, one-dimensional WENO reconstruction that we have already carried out within the coarse zone, we already have all the terms that go into eqn. (2.19). If the curl is non-zero we can also reconstruct the curl within the zone shown in Fig. 5  V at the four edges that lie interior to Fig. 5; where the superscript "FV" indicates that they were obtained from the two-dimensional finite volume WENO reconstruction. Clearly, because the multidimensional WENO reconstruction on the coarse mesh was not based on any constraint-preserving principles, we do not expect the 11 11 1 1 to be equal to 1 z R in the fine mesh zone 1 shown in Fig. 5. Similarly, we do not expect Now think back to the important insight that we derived from Table I. We realized that even if the curl constraints are not exactly satisfied, they are indeed quite closely satisfied. This means that we should be able to make a small redistribution of the values ; ; x y xc yc z We can add an additional equation that allows us to parametrize the solution in terms of the variable "α " as follows:- Now notice that the four equations in eqns. (6.2) and (6.3) can be solved in terms of the parameter "α " as follows:- In the above equations, the linear dependence on the parameter "α " has been made explicit. But how do we set the parameter "α "? We set it by demanding that the set of variables The final value that we get for the parameter "α " after the least squares minimization is given by Substituting the above value for the parameter "α " in eqn. (6.4) gives us the optimal solution that satisfies the discrete circulation conditions exactly while remaining as close as possible to the higher order values from the multidimensional, finite volume WENO reconstruction. This completes the process of showing how the curl-preserving prolongation can be carried out in two dimensions with a refinement ratio of 2.
There is still one very important observation that we make from the procedure we have described above. Notice that a specification of the values of the values for V at the internal edges. Of course, at second order, the description in Sub-section II.2 is much easier to implement. However, the balance tips over as we go to substantially higher orders because the description in this Sub-section simply uses a standard multidimensional, finite-volume WENO reconstruction. Notice that at higher orders, the values FV yc V will be responsive to the higher order FV reconstruction. Therefore, the parameter "α " will indeed be different at each increasing order. In other words, the touch-up procedure described here is truly responsive to the order of the FV reconstruction. Fig. 6a shows the situation where we have a refinement ratio of 3. In that case, we have 12 internal edges where we can specify the vector components using a high order, two-dimensional, finite volume reconstruction. These locations are shown by the red arrows. There are 8 independent curl conditions which can be used to obtain parametrically-defined improved solutions at the location of the blue arrows. This implies a solution strategy based on 12-8=4 free parameters which can be optimized to preserve the order property as closely as possible. In practice, we would use a constrained least squares optimization process for solving such problems. In Fig. 6a the subdivisions of the coarse zone are shown to be uniform because we wish to demonstrate an AMR prolongation step. However, if one wanted to use these subdivisions to carry out a higher order reconstruction, it is even possible to shift the magenta lines in Fig. 6a so that they coincide with some suitable choice of quadrature points. Fig. 6b shows the situation where we have a refinement ratio of 4. In that case, we have 24 internal edges where we can specify the vector components using a high order, multidimensional, finite volume reconstruction. These locations are shown by the red arrows. There are 15 independent curl conditions which can be used to obtain parametrically-defined improved solutions at the location of the blue arrows. This implies a solution strategy based on 24-15=9 free parameters which can be optimized to preserve the order property as closely as possible. Again, this is most easily done using a constrained least squares optimization process. As before, if one wanted to use these subdivisions to carry out a higher order reconstruction, it is even possible to shift the magenta lines in Fig. 6b so that they coincide with some suitable choice of quadrature points. The pattern is now easy to generalize. If we have a refinement ratio of "n", we would have ( ) 2 1 n n − internal edges and 2 1 n − independent curl conditions with the result that we would have ( ) 2 1 n − free parameters which can be used to preserve the order property as closely as possible.
Asymptotically, as n → ∞ , the ratio of the number of internal edges to the number of free parameters tends to 2; showing that there is considerable freedom in the number of free parameters even at high refinement ratios or high orders of reconstruction. This small calculation shows us that the procedure designed in this Sub-section can be extended to all orders.

IV.2) Restoring Curl-Preserving Structure in Three-Dimensions
The narrative in the previous Sub-section is indeed very useful even when carrying out curl-preserving prolongation on three-dimensional AMR mesh hierarchies. This is because we initialize edges of the second type on the fine mesh, i.e. the edges on the fine mesh that lie within the faces of the coarse mesh, before we proceed to initializing edges of the third type. However, this still requires us to consider edges of the third type on the fine mesh, i.e. the edges of the fine mesh that lie entirely within a coarse mesh zone. We consider those types of edges in this Subsection.
In Fig 7, the edge numbering is loosely patterned after Fig. 3. As in Fig. 5, the solid black lines show one zone of a coarse 3D mesh; while the dashed magenta lines demarcate the eight fine mesh zones that are obtained by sub-dividing the coarse mesh zone with a refinement ratio of 2. As in Fig 5, at each edge of the coarse mesh, a one-dimensional higher order WENO reconstruction is applied to the vector field component in the direction of that edge. This is used to obtain two higher order interpolated values along the two coincident edges of the fine mesh, which gives us the type one initialization. Within each coarse mesh face, we use the construction from the previous Sub-section to obtain unique edge values for the fine 3D mesh, which gives us type two initialization (see equations (6.2) to (6.6)). For example, in the xz-face nearest to us, we obtain . These edge-centered values remain the same across adjoining zones of the coarse mesh, which gives us reason to think that a globally curl-preserving touch-up strategy can be invented for the fine mesh in three dimensions.
We now focus on the procedure for initializing the type three edges of the fine mesh. We begin by obtaining zone-averaged values for all three components of the vector field on the coarse mesh using eqn. (3.29). For each of the three zone-averaged components, we carry out a high order, three-dimensional finite-volume WENO reconstruction. This reconstruction allows us to obtain the high order accurate red values along the newly-introduced edges of the fine mesh. In practical terms, this allows us to specify 1  Fig. 7; where the superscript "FV" indicates that they were obtained from the multidimensional finite volume WENO reconstruction. These values are shown as red arrows and labeling in Fig. 7. These vector field components are order-preserving on the fine mesh, but they are certainly not curl-constraint preserving. The task at hand is to replace them with their nearest approximations which are indeed curl-constraint preserving. The curl-constraint preserving components of the vector field that are newly-introduced on the fine mesh are shown with blue arrows and labeling in Fig. 7. They are denoted by 1 Fig. 7. The insight that we have gained from Table I give us assurance that a small touch up of 1  Fig. 7).
The faces of the fine mesh that are internal to the coarse mesh zone are shown in Fig. 7. We can see that there are 12 such fine mesh faces. However, the circulation is not independent in these 12 fine mesh faces. Indeed, using arguments that are similar to those developed with eqns. V and we label that parameter " β ". To find the optimal value of " β ", we will use a least squares minimization to nudge these values as close as possible to the order preserving (but unconstrained) values given by 1 Since the mathematical steps are entirely analogous to the ones that are explicitly shown in the previous Sub-section, we will not go through the steps here. We will just provide the final answer that can be implemented in code. The parameter " β " which optimizes the solution is given by 11 This optimal parameter can then be used to obtain the optimal values for  x x The round brackets in eqns. (6.7) and (6.8) should not be thought of as matrices. This completes the process of showing how the curl-preserving prolongation can be carried out in three dimensions with a refinement ratio of 2.
When we have a refinement ratio of 3, we have 36 internal edges where we can specify the vector components using a high order, multidimensional, finite volume reconstruction. There are 28 independent curl conditions which can be used to obtain parametrically-defined improved solutions at the location of the blue arrows. This implies a solution strategy based on 36-28=8 free parameters which can be optimized to preserve the order property as closely as possible. In practice, we would use a constrained least squares optimization process for solving such problems. For a refinement ratio of "n", we would have free parameters which can be optimized to preserve the order property as closely as possible. Asymptotically, as n → ∞ , the ratio of internal edges to the free parameters tends to 3; which shows that we continue to have considerable freedom in the number of free parameters even at high refinement ratios or high orders of reconstruction. The subdivisions of the coarse zone can be uniform if we are carrying out an AMR prolongation step.
However, if one wanted to use these subdivisions to carry out a higher order reconstruction, it is even possible to make the subdivisions non-uniform so that they coincide with some suitable choice of quadrature points.

V) Results: von Neumann Stability Analysis of Curl Constraint-Preserving DG Schemes
A von Neumann stability analysis for classical, volume-centered, DG schemes has been done (Liu et al. [47], Zhang and Shu [62]). In that analysis, the authors focused on the advection equation with a constant velocity. In a classical DG scheme, the primal variables are zone-centered and the objective is to satisfy a telescoping density-conservation constraint. (In other words, mass, momentum and energy densities are conserved in any subset of zones because of a telescoping application of mass, momentum and energy fluxes.) A similar von Neumann stability analysis for the induction equation in MHD has also been carried out by Balsara and Käppeli [20]. The induction equation evolves a vector field in divergence-free fashion, and with the simplification of a constant velocity in two dimensions, it also reduces to an advection of the two components of the vector field. In a divergence-preserving, face-centered, DG-like scheme for the induction equation, the primal variables are face-centered and the objective is to satisfy a telescoping divergence-preserving constraint. This choice of collocation also holds true for any divergenceconstraint preserving PDE like Maxwell's equations (Balsara and Käppeli [24]). It is now easy to see that for the curl-constraint preserving, edge-centered, DG schemes we have the objective to satisfy a telescoping curl-preserving constraint. As mentioned before, with a constant velocity, eqns. (2.23) show us that the model PDE again reduces to the advection of two curl-free vector components. In all such stability analyses it is traditional to simplify the spatial and temporal parts of the problem by using multi-stage Runge-Kutta timestepping. Therefore, we will use the SSP-RK timestepping schemes from (Shu and Osher [56], [57], Shu [58], Spiteri and Ruuth [54], [55], Gottlieb et al. [43]). The temporal order in our von Neumann stability analyses will always be matched to the spatial order of accuracy of the DG scheme.
We saw in Sub-sections II.5 and III.5 that the ingredients of a successful curl constraintpreserving scheme consist of :-1) a higher order curl-free reconstruction, 2) a multidimensional Riemann solver and 3) a suitable high order timestepping strategy. It is possible to use these three building blocks to design Finite Volume, PNPM and DG schemes of higher order Godunov type that preserve the global constraints. The full description of such a plan requires indeed a separate paper and such a paper is under construction (Balsara and Käppeli [28]). It would take us too much out afield to provide all details here, but we mention that such a scheme, as designed in twodimensions, uses the moments in the edges of the mesh as the primal variables for the curlpreserving vector field. From those primal variables, the reconstruction strategy described in Section II is used to reconstruct the entire vector field. (The von Neumann stability for fourth order accurate, curl-preserving DG schemes has so far proved analytically intractable.) The application of the multidimensional Riemann solver, along with the use of a suitably higher order SSP-RK timestepping scheme then completes the update strategy. In Balsara and Käppeli [28] we present a von Neumann analysis of such a globally curl-free DG scheme, along with its FV and PNPM variants. In this section, as part of our results, we show some of the output from the von Neumann stability analysis applied to a globally curl-free DG scheme. This is done in the spirit of illustrating the value of curl-preserving reconstruction in numerical scheme design.
In such a von Neumann stability analysis, one posits a curl-free vector field in 2D with wave vector ( ) What matters is The amplification factor and phase depend on:-1) the relative angle between the wave vector and the velocity and 2) the ratio of the wavelength to the mesh size. In our von Neumann stability analysis we used a 2D Cartesian mesh with square zones. We consider situations where the velocity vector makes angles of 0 o , 15 o , 30 o and 45 o relative to the xdirection of the mesh. This gives us a sufficiently interesting range of velocity directions. For each of those velocity directions we allow the wave vector to sweep over all possible angles between the direction of the velocity and the direction of the wave vector. DG schemes usually display very good resolving capabilities, so we display the amplitude and phase information when the wavelength spans 5 zones, 10 zones and 15 zones. We do this for P=1 and P=2 DG-like schemes, with the result that we consider second and third order accurate DG-like schemes. Despite our having adopted the "DG" nomenclature, please note that the globally curl-free DG-like schemes discussed here are very different in construction from the classical DG schemes. Fig. 9 shows the wave propagation characteristics for globally curl-free P=1 DG schemes, which are second order accurate. In these schemes, each edge has two time-evolutionary modes for the vector field, but the modes are connected because of the curl-free condition. Figs. 9a to 9d show one minus the absolute value of the amplification factor when the velocity vector makes angles of 0 o , 15 o , 30 o and 45 o relative to the x-direction of the mesh. Figs. 9e to 9h show the phase error, again for the same angles. The 2D wave vector can make any angle relative to the 2D direction of velocity propagation, with the result that the amplitude and phase information are shown w.r.t. the angle made between the velocity direction and the direction of the wave vector. In each plot, the blue curve refers to waves that span 5 cells per wavelength; the green curve refers to waves that span 10 cells per wavelength; the red curve refers to waves that span 15 waves per wavelength. Notice that the cases where the velocity vector makes angles of 0 o and 45 o relative to the x-direction of the mesh indeed show symmetrical wave propagation characteristics, as expected. When the velocity vector makes angles of 15 o and 30 o relative to the x-direction of the mesh, there is no symmetry between the velocity direction, the mesh direction and the direction of the wave vector, with the result that we don't expect to see symmetrical plots, and indeed we don't. We do, however, observe that when the waves span 10 cells per wavelength and 15 cells per wavelength the wave propagation becomes very close to isotropic and quite free of dissipation. This is a good sign that even our second order DG scheme shows rather isotropic wave propagation with increasing wave length. In all instances, Figs. 9a to 9d show us that one minus the amplification factor is always positive or zero, indicating that the globally curl-free, second order, DG scheme is indeed stable. Fig. 10 shows the wave propagation characteristics for globally curl-free P=2 DG schemes, which are third order accurate. In these schemes, each edge has three time-evolutionary modes for the vector field, but the modes are connected because of the curl-free condition. The first four panels in Fig. 10 show one minus the amplification factor, and the next four panels in Fig. 10 show the phase error. Comparing these results to their analogues in Fig. 9 we see that there has been an order of magnitude improvement in the amplification factor as well as the phase. This shows the benefit of resorting to a higher order DG scheme. We also see that the wave propagation in Fig.  10 is much more isotropic compared to Fig. 9. This shows that in multidimensions, higher order DG schemes provide not just vastly improved accuracy but also significantly improved isotropy of wave propagation.
Taken together, the results of this Section show that our curl constraint-preserving reconstruction, coupled with the multidimensional Riemann solver, provides a successful framework for the design of a very proficient class of mimetic DG-like schemes for involutionconstrained PDEs. Most importantly, the results presented point to a class of DG-like mimetic schemes for involution-constrained PDEs that have superior amplitude preservation and phase accuracy even in multiple dimensions.

VI) Numerical Results for a Model Problem -Demonstration of Order Property
It behooves us to design and display a model problem where one can palpably witness the value of an exactly curl-preserving scheme. Moreover, since we have presented entire classes of such schemes with increasing order of accuracy, we would like to demonstrate the value of high order of accuracy. To that end, we present a test problem for a simple model PDE system, where the curl-free evolution of the vector field is of crucial importance.
Let us begin our discussion by considering the currently-available alternatives. Of course, a GLM-style cleaning procedure has been developed in Dumbser et al. [39], but it requires increasing signal speed of the cleaning equations and also adds many more vector fields than are originally necessary. In Dumbser et al. [40] and Boscheri et al. [31] a new exactly curl-free semiimplicit scheme was presented using a vertex-based staggered mesh, but it is limited to second order of accuracy and requires frequent interpolation of the velocity field and of the curl-free vector field J to different staggered locations on the mesh, which makes it more difficult to extend to adaptive mesh refinement (AMR) techniques.
By contrast, the present formulation preserves the same control volume for the fluid variables as well as the curl-constraint preserving vector field, making it suitable for an eventual future extension to AMR. The availability of higher order curl-preserving formulations also allows us to show another interesting facet that has gone unappreciated in the literature. (However, see Balsara [5], [8] where this same consideration has been demonstrated for divergence-preserving schemes.) It turns out that in certain important limits, the curl-preserving vector field (quite like its divergence-preserving counterpart) satisfies an energy principle. The quadratic energy of the vector field should remain unchanged in time. A good scheme should preserve this quadratic energy as much as possible. We show in this section that our increasingly accurate curl-preserving schemes preserve the quadratic energy with increasing precision.
In this Section we illustrate the capabilities of our new high order accurate numerical method with the help of the toy system introduced in Dumbser et al. [39], [40]. The governing PDE system reads where the Einstein convention implying summation over repeated indexes is adopted. The system of equations (6.1)-(6.3) describes the evolution of a scalar quantity ρ and two vector fields v and J, that in a fluid dynamic context could be interpreted as density, velocity and a kind of thermal impulse, respectively, while p represents the equivalent of a pressure, see Dumbser et al. [36]. As shown in Dumbser et al. [40] the system satisfies the extra energy conservation law In this paper the model system is closed by the simple linear relation where γ is a given constant, as well as c0 in eqn. (6.2). The system (6.1)-(6.3) with (6.4) falls into the larger class of symmetric hyperbolic and thermodynamically compatible (SHTC) systems studied by Godunov and Romenski, see Godunov [41] and Romenski [51] and references therein.
for all times if the field J was curl-free at the initial time.
It is therefore crucial to satisfy this constraint even at the discrete level, that is if 0 mk C = at the initial time it must remain zero for all times.

VI.1) Model Problem: A Stationary curl-free solution
To verify that the novel curl constraint-preserving scheme is able to fulfill over time the involution constraint in eqn. (6.6) we propose to solve the following test problem, which is an exact stationary and smooth solution of the model system (6.1)-(6.3). Let the computational domain be the square Ω=[-5;5] 2 and let the generic radial coordinate r satisfy r 2 =x 2 +y 2 . The quantity J is defined as the gradient of a scalar potential φ, so that the initial condition is ensured to be curl-free. The potential is given by with the parameters A, R0 and σ. Then, the initial condition for the radial component In Figure 11 we show the numerical results for the vector component Jx obtained on a mesh with 50 x 50 cells by running four different schemes: the semi-implicit second order accurate staggered curl-free (SCF) scheme proposed in (Boscheri et al. [31]) and the second, third and fourth order accurate edge centered curl-preserving (ECCP) reconstruction methods developed in Section II of this work. The final time of the simulation is very large, i.e. tend=100, in order to show the behavior and the stability of the scheme for very long time computations. The less dissipative behavior achieved by the high order order reconstructions is clearly visible. We also notice that the SCF scheme is the most dissipative of the schemes shown in Fig. 11 because of the copious interpolation of the velocity field and the curl-free vector field J to different staggered locations on the mesh. We also mentioned that the quadratic energy of the vector field should, in principle, be preserved for this physical problem. All numerical schemes fall short of this ideal goal. In Fig.  12 we plot out the mesh-integrated quadratic energy, simply evaluated as 2 2 2 x y J J J = + as a function of time. We see that the higher order schemes preserve the quadratic energy much better than the lower order schemes. This is because the numerical viscosity is significantly reduced when a high order reconstruction technique is adopted. Finally, we would also like to point out that without a curl-preserving scheme the solution is spoiled after short times and the equilibrium is violated very soon leading to catastrophically unphysical results. At late times, the result of not treating the curl-preserving aspect of the PDE is indeed a code blowup.
A numerical convergence study for this test problem is carried out by solving the model PDE system with the stationary initial condition discussed above until a final time of t=10 using a sequence of successively refined meshes. The results for our novel curl-free WENO schemes presented in this paper are shown in Table II for nominal orders of accuracy from two to four. In Fig. 13 we finally show the temporal evolution of the curl error of second, third and fourth order curl-free WENO schemes on a uniform Cartesian mesh composed of 32 x 32 elements. It can be clearly seen that in all cases and for all times, the curl errors remain at the level of machine precision, as expected. Taken together, the results of this Section show that our curl constraint-preserving reconstruction, coupled with the multidimensional Riemann solver, provides a successful framework for the design of mimetic finite volume schemes of increasing order of accuracy. Moreover, these schemes preserve the curl-constraint during very long time integrations. This constraint-preservation also contributes significantly to the enhanced stability of the scheme. The model problem that we provide here is also quite novel and it enables us to precisely document that the methods presented here do indeed meet their designed order of accuracy. The utility of mimetic schemes with high accuracy is also emphasized by the fact that additional quadratic energy terms are also preserved with superlative precision as one goes to higher order. As a result, we have presented high order mimetic finite volume-type schemes which have long time stability and excellent preservation of quadratic energy.

VII) Numerical Results Curl Constraint-Preserving Prolongation on AMR Meshes
In Section IV we presented a curl constraint-preserving prolongation strategy that is useful on AMR mesh hierarchies. The method was based on using eqn. (4.1) to set up a curl-free vector field on a coarse mesh and prolonging it to a fine mesh. As shown in Table I, a naïve, finitevolume-based prolongation will not preserve the curl-free structure of the vector field on the finer mesh. However, the rapidly diminishing discrete circulation on progressively refined meshes gave us the critical insight that a small touch-up procedure can be invented to mildly redistribute the vector components within the edges of the fine mesh so that we can make the solution globally curl-free on the fine mesh. Such a strategy was documented for prolongation on two-dimensional AMR meshes in Sub-Section IV.1. It was extended to prolongation on three-dimensional AMR meshes in Sub-Section IV.2. We also note that the two-dimensional strategy in Sub-Section IV.1 plays an essential role in the three-dimensional strategy from Sub-Section IV.2. Therefore, documenting a fully working curl-free prolongation in three dimensions is also an ipso facto proof that the two-dimensional curl-preserving prolongation strategy also works.
We now redo the same test that was documented in Table I, however, this time we use the methods designed in Section IV. A refinement ratio of 2 is used here. The extension and generalization of this method to all refinement ratios is deferred to a subsequent paper (Balsara and Subramanyan [29]). Table III repeats Table I, but with the additional application of the methods designed in Section IV. The discrete circulation is also shown, although it is close to machine accuracy in all instances. This shows that the methods from Section IV work. We also see that the methods all achieve their design accuracies. If anything, we see that with the touch-up procedure, the accuracies are dramatically improved for all orders in Table III as compared to  Table I. This shows that the optimization process in eqns. (6.6) and (6.7), when applied to a constraint-preserving vector field, has the extremely salutary effect of pushing the result towards increasing accuracy! While this is apparent in the third and fourth order results, it is most compellingly demonstrated in the second order results. Please compare the second order results from Table III to the second order results from Table I. We see that the imposition of curlconstraints has improved the solution by an entire order of magnitude, turning an initially second order result into a third order accurate result! Similar salutary trends are seen to a slightly lesser degree by comparing the fourth order results in Table III to the fourth order results in Table I. Table III shows the result of using the touch-up methods from Section IV along with standard finite volume WENO reconstruction at second, third and fourth order to prolong a curl-free vector field from a coarse mesh to a finer mesh with a refinement ratio of two.
Only the x-component of the curl-free vector field is shown because the other components show a similar trend. We show the errors in the prolongation of the vector field. We also show the maximum circulation of the vector field on the faces of the fine mesh with the intention of documenting that it is close to machine precision. We have also made the case that the constructive process in Section IV can also be used as an alternative strategy for curl-preserving reconstruction. If that claim is to be properly demonstrated, we should be able to take the second and third order reconstruction strategies from Sub-sections III.2 and III.3 respectively and use them to carry out curl-preserving prolongation on AMR meshes. In other words, the availability of an analytically exact curl-preserving reconstruction from Sub-sections III.2 and III.3 can also be used as an alternative strategy for curlpreserving prolongation on AMR meshes. If the touch-up-based procedure from Section IV is good, then it should yield results that are competitive with the analytically exact curl-preserving results from Sub-sections III.2 and III.3. Table IV shows such results where Sub-sections III.2 and III.3 were used to prolong the constrained vector field with second and third order accuracy. Comparing Tables III and IV, we see that the results are indeed entirely competitive! From Table  IV we see that the second order results are effectively third order accurate. This is only true for refinement ratios of 2 because there is a serendipitous cancellation of third order terms when a refinement ratio of 2 is used. It would not be true for other refinement ratios. Table IV shows the result of using the analytically exact curl-preserving methods from Subsections III.2 and III.3 at second and third order to prolong a curl-free vector field from a coarse mesh to a finer mesh with a refinement ratio of two. Only the x-component of the curl-free vector field is shown because the other components show a similar trend. We show the errors in the prolongation of the vector field. We also show the maximum circulation of the vector field on the faces of the fine mesh with the intention of documenting that it is close to machine precision. Taken together, the results from this Section show that we have indeed arrived at two equivalent strategies for curl constraint-preserving reconstruction. The first one is analytical and explored in Sections II and III. The second one builds on traditional finite volume reconstruction and uses a small touch-up procedure to restore the constraints. The latter strategy may be deemed more accessible than the former by most practitioners, though the former strategy has a conceptual elegance of its own. The analytical approach was also critically useful for the von Neumann stability analysis presented in Section V, because the approach from Section IV would never have allowed us to carry out such a von Neumann stability analysis.

VIII) Conclusions
Structure-preserving PDEs are becoming increasingly important for several applications in physics and engineering. Of particular interest in this paper are PDEs that preserve the curl of a vector field. Examples of such PDEs include hyperelasticity, compressible multiphase flows with and without surface tension, first order reductions of the Einstein field equations as well as the novel first order hyperbolic reformulation of the Schrödinger equation, to name a few examples. We have designed methods in this paper for increasingly high order numerical treatment of such PDEs. The essential building blocks for such methods are shown to be a curl constraint-preserving, high order accurate reconstruction strategy and the use of multidimensional Riemann solvers that are needed for properly upwinded constraint-preserving time update.
Sections II and III show how curl-preserving reconstruction is carried out in two and three dimensions. In those sections we also demonstrate that when the reconstruction is combined with multidimensional Riemann solvers, we get numerical schemes that are multidimensionally stable.
The solution of such PDEs also needs to be carried out on adaptive mesh refinement (AMR) hierarchies. The most challenging problem in developing AMR for structure-preserving PDEs consists of making a constraint-preserving prolongation of a coarse mesh vector field to a fine mesh. To that end, we start with a naïve, finite-volume-based prolongation which can be carried out at any order. To that, we add a computationally inexpensive touch-up procedure to mildly redistribute the vector components within the edges of the fine mesh so that we can make the solution globally curl-free on the fine mesh. This invention is documented in Section IV.
Section V presents some of the most impressive results of a von Neumann stability analysis of globally structure preserving DG-like schemes in multiple space dimensions. The results in Section V point to a class of DG-like mimetic schemes for involution-constrained PDEs that have superior amplitude preservation and phase accuracy even in multiple dimensions.
In Section VI a test problem is constructed that produces steady-state analytically exact solutions. In that Section we show that the mimetic finite volume schemes that use our methods indeed preserve order of accuracy while simultaneously satisfying the curl-free constraint. In some limits, these schemes also preserve the quadratic energy on the mesh. The utility of our mimetic schemes with high accuracy is also illustrated by the fact that additional quadratic energy terms are preserved with superlative precision as one goes to higher order. As a result, we have presented high order accurate mimetic finite volume-type schemes which have long time stability and excellent preservation of quadratic energy.
Section VII shows curl constraint-preserving prolongation of vector fields on AMR hierarchies using the methods developed in Section IV. The extreme efficiency of our curlpreserving prolongation method stems from the fact that it draws on well-worn and well-weathered higher order finite volume reconstruction. To that we simply add a touch-up procedure to restore the constraints that are satisfied by the vector field. An unanticipated, but happy, consequence is that the touch-up procedure can actually improve the accuracy of the prolongation! In summary, in this paper we have shown that curl-constraint preserving reconstruction provides a rich set of insights for mimetic scheme design when dealing with PDEs that have a curlpreserving involution.
matched. We also see that there is sufficient flexibility with setting xx c , zz a to match the terms that are linear in "x" and "z" respectively to any desired modal values. Evaluating ( ) z As before, we look for solutions with xzz yzz b a = . We see that the terms that are linear in "z" will indeed match with the corresponding terms in eqn. (3.9). The above two sentences ensure that the values of the discrete circulations at the top and bottom z-faces of the cube in Fig. 3 will be matched. We also see that there is sufficient flexibility with setting xx b , yy a to match the terms that are linear in "x" and "y" respectively to any desired modal values.
In this paragraph we provide explicit expressions for the curl of the third order accurate vector field that is catalogued in eqns. (3.18), (3.19) and (3.20).