Efficient WENO‐Based Prolongation Strategies for Divergence‐Preserving Vector Fields Dinshaw S. Balsara1,2 · Saurav Samantaray1 · Sethupathy Subramanian1

Adaptive mesh refinement (AMR) is the art of solving PDEs on a mesh hierarchy with increasing mesh refinement at each level of the hierarchy. Accurate treatment on AMR hierarchies requires accurate prolongation of the solution from a coarse mesh to a newly defined finer mesh. For scalar variables, suitably high-order finite volume WENO methods can carry out such a prolongation. However, classes of PDEs, such as computational electrodynamics (CED) and magnetohydrodynamics (MHD), require that vector fields preserve a divergence constraint. The primal variables in such schemes consist of normal components of the vector field that are collocated at the faces of the mesh. As a result, the reconstruction and prolongation strategies for divergence constraint-preserving vector fields are necessarily more intricate. In this paper we present a fourth-order divergence constraint-preserving prolongation strategy that is analytically exact. Extension to higher orders using analytically exact methods is very challenging. To overcome that challenge, a novel WENO-like reconstruction strategy is invented that matches the moments of the vector field in the faces, where the vector field components are collocated. This approach is almost divergence constraint-preserving, therefore, we call it WENO-ADP. To make it exactly divergence constraint-preserving, a touch-up procedure is developed that is based on a constrained least squares (CLSQ) method for restoring the divergence constraint up to machine accuracy. With the touch-up, it is called WENO-ADPT. It is shown that refinement ratios of two and higher can be accommodated. An item of broader interest in this work is that we have also been able to invent very efficient finite volume WENO methods, where the coefficients are very easily obtained and the multidimensional smoothness indicators can be expressed as perfect squares. We demonstrate that the divergence constraintpreserving strategy works at several high orders for divergence-free vector fields as well as vector fields, where the divergence of the vector field has to match a charge density and its higher moments. We also show that our methods overcome the late time instability that has been known to plague adaptive computations in CED.


Introduction
Higher order Godunov methods for the stable and robust simulation of hyperbolic systems have been presented in the literature (Van Leer [40], Colella and Woodward [18], and many more). Soon after the vigorous development of such methods, one saw the emergence of adaptive mesh refinement (AMR) for these methods at second order of accuracy (Berger and Oliger [13], Berger and Colella [12]). The emergence of weighted essentially nonoscillatory (WENO) methods (Harten et al. [26], Shu and Osher [36,37], Liu et al. [31], Jiang and Shu [27], Balsara and Shu [9], and many more) showed that higher order accurate methods could indeed be designed for the treatment of fluid flow in particular and hyperbolic systems in general. These methods have subsequently led to AMR methods for fluids that operate with higher order of accuracy (McCorquodale and Colella [34], Dumbser et al. [23]). Especially when it comes to higher order methods, an important step in the implementation of an AMR solution strategy consists of the prolongation of the solution from one mesh in the refinement hierarchy to the next finer mesh in the refinement hierarchy. Prolongation is the act of accurately transferring the solution from a coarse mesh to a newly defined fine mesh. Without such an accurate prolongation, the overall accuracy is impossible to maintain in an AMR simulation.
In parallel to the above-mentioned development, there was an effort to numerically simulate PDEs with involution constraints. These are PDEs, where the structure of the PDE is such that an additional constraint is automatically preserved in the time-evolution of the PDE. Maxwell's equations, which form the basis for computational electrodynamics (CED), are prototypical of a PDE system that maintains an involution constraint. Maxwell's equations, which are based on a curl-type update, ensure that the magnetic induction remains divergence-free forever, whereas the divergence of the electric displacement is always equal to the local charge density. The Yee scheme (Yee [41]), which later led to the development of the finite difference time domain (FDTD) method for CED (Taflove and Brodwin [38], Taflove and Hagness [39]), was the first successful exemplar of a scheme that maintained the discrete divergence constraints for the electric and magnetic fields up to machine precision on a finite difference mesh. Much of the success of FDTD derives from the fact that the globally divergence-preserving constraints that are inherent in Maxwell's equations are maintained at the discrete level in FDTD. Figure 1a shows the staggering of vector field variables in FDTD which provides a direct interpretation of the two curl-type equations given by Faraday's Law and the generalized Ampere's Law, and a natural satisfaction of the constraint equations given by Gauss's Laws for electric and magnetic charge. Notice that the E (electric) and H (magnetic) vector fields in Fig. 1a form two staggered control volumes, giving rise to an automatic preservation of the divergence constraints in FDTD. However, the need for two staggered control volumes, instead of a single control volume, can act as an impediment to the development of the most convenient types of AMR methods.
The equations of magnetohydrodynamics (MHD) can be viewed as another example of an involution constrained system. Since the MHD equations rely on Faraday's law from the original Maxwell equation set (along with a constitutive relation for the electric field), the magnetic induction remains divergence-free throughout its time-evolution. It has also been shown by Brackbill and Barnes [16] and Brackbill [15] that the presence of a divergence in a numerical code can result in fictitious forces on the magnetized plasma. Drawing from the Yee scheme, Brecht et al. [17], Evans and Hawley [24], DeVore [22] designed globally divergence-free methods for MHD that were not based on higher order Godunov technology. However, part of the MHD equation set mirrors the Euler equations, for which higher   Fig. 1 a Collocation of electric fields E, and magnetic fields H, for the FDTD scheme for CED. Notice the need for two staggered control volumes which enable divergence-constraint-preserving curl-type update equations for both the vector fields. b Which applies to the MHD equations, shows the collocation of the face-centered magnetic induction, B, which is updated using the multidimensionally upwinded electric fields, E, at the edges of the control volume. The fluid variables in the MHD system have zonecentered collocation. c Which applies to the FVTD-based update of Maxwell's equations, shows the facecentered electric displacement and magnetic induction, i.e., the vector fields D and B, respectively. These are updated in a way that preserves the global divergence using a discrete circulation of the edge-centered magnetic and electric fields, i.e., H and E, respectively. Notice that all variables in b and c are defined on the same control volume Balsara [1]. For constrained vector fields, prolongation is the act of accurately transferring the vector field from a coarse mesh to a newly defined fine mesh while respecting the constraints on the finer mesh. The second-order accurate, divergence-free prolongation of vector fields was found to be of crucial importance in the design of such methods. In fact, the divergencefree reconstruction of the magnetic vector field within each zone was the crucial bottleneck in the prolongation problem, as well as the crucial bottleneck in the development of secondorder accurate methods for AMR-MHD in general. It was soon realized (Balsara [2]) that the divergence-free reconstruction of the magnetic field could play an important role in higher order scheme design for MHD. WENO methods played an important role in the development of these ideas, because the facial moments of the magnetic field had to be reconstructed using higher order WENO philosophy. The emergence of multidimensional Riemann solvers (Balsara [3,4]) was another crucial building block for numerical MHD. These advances subsequently led to globally constraint-preserving methods for CED that were based on higher order Godunov philosophy (Balsara et al. [6]). Because those methods drew on finite volume-type methods, they were called finite volume time-domain (FVTD) methods. In those papers, the globally constraint preserving reconstruction methods were extended so that the divergence of the electric displacement vector field could match all the higher order moments of the charge density. In other words, the reconstruction strategy had to be upgraded so that it was divergence-preserving rather than divergence-free. Figure 1c shows the collocation of variables in a globally constraint-preserving higher order FVTD scheme for CED. Comparison with Fig. 1a shows that these more modern methods are defined on a single control volume-opening the door to very easy implementation of AMR. Progress in higher order globally divergence constraint-preserving MHD and CED schemes has been swift since then. Globally divergence constraint-preserving discontinuous Galerkin-like (DG) methods for MHD and CED systems were designed and analyzed (Balsara and Käppeli [7]) and their implementation was documented in Hazra et al. [25], Balsara et al. [8].
The goal of this paper is to design methods for divergence-preserving prolongation of vector fields at high order. Furthermore, we are only interested in the three-dimensional case, because the two-dimensional problem is not of much interest in practical AMR applications. By default, unless it is specified, we will consider refinement ratios of two. However, the methods are general and, in the later sections, we will show that they can be used for refinement ratios that are larger than two. At second order, the problem was solved in Balsara [1] who presented a polynomial-based reconstruction strategy that could be used for prolongation. We present a very brief synopsis of that strategy so that the reader can appreciate the options available to us as we try to push towards higher order. Figure 2 shows a typical situation, where a fine mesh abuts a coarse mesh. If the coarse mesh has to be refined, we require that all the information about the four vector field components from the adjoining four fine mesh faces should be retained in the vector field that is reconstructed in the abutting coarse mesh zone. If we retain just the four components that are present in the four fine mesh faces, then it means that each coarse mesh zone, where prolongation is to be carried out should be able to accommodate up to four pieces of information at each of its six faces. Using this information, one has to find three higher order volume-filling polynomials for the three vector field components in the coarse zone that is about to be refined. The polynomials should be such that they can match up to four pieces of information at the six faces of the coarse zone; they should do so while satisfying all the divergence-free constraints. This creates quite a mathematical puzzle, but at second order it was solved in Balsara [1]. (Subsequent work in Balsara et al. [11] has shown that the divergence-preserving constraints can also be accommodated.) Once the three polynomials are found, the vector field is analytically specified and it can be prolonged via simple areal integration and averaging to all the faces of any set of refined zones that replaces the coarse zone in question. Section 3.1 of Balsara [1] shows a very simple and easy to follow example of how this is done in two-dimensions and Sect. 4 of the same paper provides all the three-dimensional details. Now let us identify the challenges as we try to transition to the prolongation problem at higher order. We will include the consideration of situations, where a DG-like scheme might be involved. Even the simplest DG-like scheme that retains first moments of the vector fields in the faces will indeed retain a mean value plus two first moments of the vector field components in each of the four fine faces shown in Fig. 2. Consequently, the three reconstructed constraint-preserving polynomials within an abutting coarse zone that requires refinement have to be able to match at least twelve pieces of information at each of its faces. This puts considerable pressure on the search and discovery of a set of three higher order polynomials that can do that. Even if one just wishes to design a prolongation strategy for an FVTD-like algorithm at high order (say third or fourth order), it still implies that a lot of facial complexity from one coarse zone has to be transferred in an order preserving fashion to the faces of a refined set of zones. In Sect. 2 of this paper, along with its associated supplements, we present an analytically exact, fourth order accurate strategy for carrying out divergence constraint-preserving prolongation of vector fields. In other words, despite the complexity of the polynomial space that has to be explored, the task proves tractable.
Higher order WENO reconstruction of the facial components of the vector field is a first step in the process that is described in the above paragraph, and in this paper we have presented some very novel innovations on that front. In particular, Supplement A of this paper Fig. 2 Fine mesh that abuts a coarse mesh. A refinement ratio of two is shown with the result that four fine mesh faces abut a coarse mesh face in the interfacial region. The green dashed region shows a hypothetical new fine mesh that has to be initialized. Some of the zones in the green region can be initialized by straight injection from the existing fine mesh. A straightforward high order reconstruction would suffice for some of the other zones. However, the layer of zones that abut the fine mesh will have to retain all the moments that are present in the faces of the fine mesh. They need special attention and the analytically exact reconstruction presented here is designed to accommodate such cases shows that finite volume WENO reconstruction can be quickly and optimally carried out on structured meshes with closed form smoothness indicators that can be expressed as the sum of perfect squares. This innovation has never been documented in the literature before, and it will be useful even to those whose interests are purely in higher order finite volume WENO-based reconstruction. To put it in perspective, in Balsara et al. [5] we were able to show for the very first time that the smoothness indicators for higher order finite difference WENO can be written as the sum of perfect squares and can be evaluated very efficiently if cast in a basis set of Legendre polynomials of increasing order. In this paper, the analogous exercise has been presented for finite volume WENO reconstruction on Cartesian meshes for the very first time. Supplement A of this paper can, therefore, have stand-alone interest for certain WENO practitioners.
Despite the advances described in the previous two paragraphs, some of the steps in the discovery of the three polynomials that make up the divergence-constraint-preserving vector field cannot be easily extended past fourth order of accuracy. The reason is that the number of terms that have to be juggled, and the number of constraints that have to be met, prove to be too large as the order of accuracy is increased. However, it is almost certain that people will want to carry out constraint-preserving prolongation of vector fields at orders that go beyond fourth order. To do this we introduce another innovation that can indeed be extended to all orders in a relatively comfortable fashion. We describe that next.
Realize that the vector fields in the faces of the coarse mesh have to accommodate a certain number of moments. As the accuracy of the base-level scheme increases, the number of facial moments increases. Therefore, one has to invent a WENO-like reconstruction strategy that matches the values at either of the two opposing faces of each zone. This WENO-like reconstruction strategy does not have to be constraint-preserving. However, it does have the following special property: even if the divergence constraint is violated, as the order of accuracy is increased, the extent to which the constraint is violated decreases in an orderly fashion with increasing mesh refinement. In other words, the third, fourth, and fifth order WENO-like reconstruction strategies that we present in Sect. 3 will nevertheless have residuals that deviate from the exact constraint condition with second, third, and fourth orders of accuracy, respectively! For this reason, we refer to the WENO-like reconstruction strategy that we present here as an almost divergence-preserving reconstruction strategy (or WENO-ADP). The invention of this almost divergence-preserving reconstruction strategy is very good news, because it indicates that the extent to which the constraint is violated is strictly controlled, even though we shall make no special effort in Sect. 3 to preserve the constraint. The implication of this insight is that we can invent a touch-up procedure that will fully restore the constraint in each zone. The touch-up procedure does not have to do much to restore the constraint exactly. When the refinement ratio is two, we provide analytically exact expressions for restoring all divergence constraints exactly within a newly refined fine mesh zone; this is presented in Sect. 4. When the refinement ratio is greater than two, we present a constrained least squares (CLSQ) based method for restoring the divergence constraint up to machine accuracy; this is described in Sect. 5. The WENO-like almost divergence-preserving reconstruction strategy with touch-up is referred to as WENO-ADPT. Furthermore, we show that this approach should work on meshes with zones that do not have to be Cartesian but rather can have any shape.
In this paper, we do not develop higher order divergence constraint-preserving prolongation for non-Cartesian meshes. Even so, it is worthwhile stressing the forward-looking import of the insight developed here with a couple of examples. Balsara and Dumbser [5] have shown that divergence-free MHD can be carried out on tetrahedral meshes. The innovation presented here, if properly developed, will one day open up the prospect of doing high order divergence constraint-preserving AMR on all such mesh geometries. We see, therefore, that proper development of the ideas presented here will have some far-reaching consequences.
Section 2 presents a fourth order accurate divergence-preserving reconstruction strategy for vector fields that can be used for fourth order prolongation of such fields. Section 3 documents the WENO-like reconstruction strategy that matches the values at either of the two opposing faces of each zone. The reconstruction does not preserve the divergence constraint exactly, but the errors that arise from the violation of the constraint are shown to be well-controlled, i.e., it is the WENO-ADP algorithm. Section 4 presents an analytically exact strategy for restoring the divergence constraint exactly when prolongation is carried out on mesh hierarchies with a refinement ratio of two. We call this a touch-up procedure. Section 5 extends this touch-up procedure to other refinement ratios. Together, Sects. 4 and 5 present the WENO-ADPT algorithm for the divergence constraint-preserving prolongation of vector fields on adaptive meshes. WENO-ADPT can, therefore, be pronounced as "WENO-adapt". Section 6 presents implementation-related details for prolongation. Section 7 explains how an AMR timestep has to be constructed for a divergence-constrained PDE system. Section 8 presents results from prolongation. Section 9 shows how our methods provide a perfect resolution of the long-time instability that has plagued CED calculations that use mesh refinement. Section 10 presents conclusions.

Analytically Exact, Constraint-Preserving Reconstruction at Fourth Order That Is Suited for AMR Prolongation
Taking CED as our motivation, we consider the electric displacement vector field, , and the charge density, , which together satisfy the divergence constraint The ideas developed can, of course, be applied to any vector field and any density. Setting the charge density to zero yields a divergence-free reconstruction strategy. While some parts of the fourth order formulation were presented in Balsara et al. [11], that formulation was not general enough to be applicable to the fourth order accurate divergence-preserving prolongation problem for AMR. The reason is evident from Fig. 2. To retain fourth order reconstruction of a facial component, one has to retain only ten moments within a face. (For example, in the z-face one would have to retain a constant term along with x-and y-moments, plus x 2 -, y 2 -, and xy-moments, and additionally the x 3 -, y 3 -, x 2 y-, and xy 2 -moments; for a total of ten moments.) However, consider a situation, where each fine mesh face had to retain not just the mean value of the normal component of the electric displacement but also the two linear moments in the two transverse directions, as shown in Fig. 2. That would mean that one has to have a minimum of twelve moments in each face. The formulation in Balsara et al. [11] does not have that extra flexibility, whereas the formulation presented here does have such flexibility.
We now present the basics of an analytically exact constraint-preserving reconstruction at fourth order, as it should be adapted to the prolongation problem. Let the zone size be Δx , Δy , and Δz in the x-, y-, and z-directions. It is more economical to present the results for the fourth order case in the space of a reference element spanning −1∕2, 1∕2 3 . At the right and left x-faces of the reference element, the fourth order accurate facial reconstruction of the x-component of the electric displacement is given by The fourth order accurate facial reconstruction for the y-component of the electric field at the upper and lower y-faces is given by The analogous facial reconstruction for the z-component of the electric field at the top and bottom z-faces is given by We wish to find a solution in the interior of the zone in question, consistent with the facial variations in Eqs. (2), (3), and (4). Please note that the second order terms in the above equations are shown in black, the third order terms are shown in red, the fourth order terms are shown in blue, and the extra terms that are needed for retaining twelve free pieces of information at each face of the coarse mesh are shown in green. The colorization of terms will prove most helpful in tracking which terms in the higher order reconstruction enter at which particular order.
We also want the solution for the electric displacement in the interior of the zone in question to be consistent with the constraint in Eq. (1). The most economical way to initialize such a computation is to transcribe D x± 0 → D x± 0 Δx and similarly for all the other coefficients in Eq. (2). For all the coefficients in Eq. (3) we make transcriptions that are analogous to D y± Δy . Likewise, for all the coefficients in Eq. (4) we make transcriptions that are analogous to D z± 0 → D z± 0 Δz . The original zone-averaged part of the divergence-preserving condition after it is subjected to this transcription becomes and the above equation is much more tractable. Once such a zone-averaged charge density is obtained in all the zones, we can apply zone-centered, finite volume WENO reconstruction to the zone-averaged charge density to get
In situations, where a coarse zone needs refinement, and its faces do not abut a fine mesh, see Fig. 2, we will still need to reconstruct up to fourth order accurate moments within such faces of the coarse zone. A WENO strategy for obtaining the ten moments that would then be needed for retaining fourth order accuracy in Eqs. (2) to (4) is described in (Balsara and Shu [9], Balsara et al. [6]). Because the description is scattered through multiple papers, we provide a compact description of the process in Supplement A, along with fifth order extensions that are genuinely novel. However, as also seen from Fig. 2, there will be some coarse zones that abut a fine mesh and are in need of refinement. In such a situation, Fig. 3 shows that we will need to retain at least twelve pieces of information to capture all the fine mesh facial information on the abutting face of the coarse mesh. Figure 3 shows the four z = constant faces of a fine mesh. Within each of those fine faces, we have a mean value and its piecewise linear variation in the two transverse directions. Using just the information contained in Fig. 3, we show that all the moments in Eq. (4) can be fully recovered. We see that the problem reduces to the inversion of a 12 × 12 linear system of equations, and the solution is given by

3
Using the procedure from the above equation when it is needed, and using plain-vanilla finite-volume WENO reconstruction when the above equation is not needed, we can always specify all the moments in Eqs. (2) to (4). Let the x-component of the electric displacement within the unit cube be described by the following polynomial: Notice that Eq. (8) has all the terms that are needed for up to fourth order accurate reconstruction. However, it also has extra terms for ensuring that the facial moments can be matched and that the constraints can be preserved. The zone-centered terms a xxy and a xxz in the above equation cannot be satisfied by just examining the values in the faces of that zone and we will specify a WENO method that looks at the x-face values of the adjoining zones in the x-direction.   Four zone-faces in one of the z = constant planes of the fine mesh that abut the z-face of a coarse mesh. Within each fine zone-face, and relative to its local facial centroid, we have a mean value as well as x-and y-moments. This means that twelve pieces of information have to be specified on the abutting z-face of the coarse mesh if it is to fully retain all the information that is available from the fine mesh 1 3 The zone-centered terms b xyy and b yyz in the above equation cannot be satisfied by just examining the values in the faces of that zone and we will need a WENO method that looks at the y-face values of the adjoining zones in the y-direction. Let the z-component of the electric displacement within the unit cube be described by the following polynomial: As before, the zone-centered terms c xzz and c yzz in the above equation cannot be satisfied by just examining the values in the faces of that zone and we will need a WENO method that looks at the z-face values of the adjoining zones in the z-direction. We see that Eqs. (8), (9), and (10) have a large number of coefficients, because those three polynomials have not just to match the facial moments at opposing faces but they also have to simultaneously satisfy all the constraints arising from Eq. (1). Supplement B provides all the details for satisfying these coefficients. Supplement B is written in a style that facilitates computer implementation.
Once the coefficients in Eqs. (8), (9), and (10) are defined, we can require that they satisfy the divergence constraint in Eq. (1) at all orders in the polynomial expansion. The need to balance all the terms in the constraint equation in a nice and symmetrical way accounts for many of the fourth and fifth order polynomial terms in Eqs. (8), (9), and (10). Please note that the black, red, and blue terms in Eqs. (8), (9), and (10) correspond to terms that are needed for second, third, and fourth order accuracy. The magenta and orange terms are terms that we have to borrow from fifth and sixth order just to ensure that all the constraints up to fourth order are fully satisfied. This satisfaction also requires us to ensure that the reconstructed polynomial minimizes the quadratic energy in the vector field. The green terms in Eqs. Recall that at the beginning of this section we divided all the coefficients of Eqs. (2), (3), and (4) by Δx , Δy, and Δz , respectively, before embarking on the constraint-preserving reconstruction described in the above paragraphs. We now undo that process so that all the coefficients in Eq. (8) are multiplied by Δx ; so that we have a 0 → a 0 Δx for example. Similarly, all the coefficients of Eq. (9) are multiplied by Δy ; so that we have b 0 → b 0 Δy for example. Likewise, all the coefficients of Eq. (10) are multiplied by Δz ; so that we have c 0 → c 0 Δz for example.

WENO-Based Almost Divergence-Preserving Reconstruction-WENO-ADP
Let us consider the simplest situation that prevails at second order as a motivating example. Say that the left and right x-faces of a coarse zone have piecewise-linear transverse moments in the y-and z-directions so that we can write them as D x− (y, z) and D x+ (y, z) in the above equation represent these second-order accurate facial variations in the left and right faces, respectively. Realize, therefore, that we are starting with the six variables, D x− 0 , D x− y , D x− z , D x+ 0 , D x+ y , and D x+ z , from Eq. (11). We ask the question: is there some polynomial that gives us D x (x, y, z) in the full three-dimensional zone such that the polynomial values match the facial variations from Eq. (11)? We realize that if we do moment-by-moment linear reconstruction in the x-direction we will get our desired answer which is We see that Eq. (12) is a full volumetric reconstruction that matches the facial variations from Eq. (11) in the left and right faces of the zone being considered. As a penalty for trying to match the entire variation in the faces, we have to retain two extra terms that vary as xy and xz , but even those terms are fewer in number than the number of terms we would have to retain if we wanted a full divergence-preserving reconstruction. Consequently, we suggest that Eq. (12) might be a good starting point for the prolongation problem at second order, if we can subsequently show that the divergence remains bounded and consistent with an order property. With this insight in hand, we now proceed to the construction of third, fourth, and fifth order WENO-ADP approaches. For the second-order case, we did not have multiple stencils to choose from. We will see that at higher order we will indeed have multiple stencils to choose from, making it valuable to adopt a WENO-like philosophy. In the next three subsections, we detail the WENO-ADP for third, fourth, and fifth orders.

r = 3 WENO-ADP
Say that the left and right x-faces of a coarse zone have up to piecewise-quadratic transverse moments in the y-and z-directions so that we can write them as We now realize that to have all the volumetric moments up to third order, we only have to do moment-by-moment linear reconstruction in the x-direction for the linear and quadratic terms in Eq. (13). As a penalty for trying to match the entire variation in the faces, we have to retain three extra terms that vary as xy 2 , xz 2 , and xyz . However, notice that we will now need to reconstruct the constant, x-dependent and x 2 -dependent variation in the x-direction. These three pieces of information cannot be obtained exclusively from D x− 0 and D x+ 0 . We will have to look to the x-face that is one zone rightward of the right face of the zone of interest to obtain D x2+ 0 . Similarly, we will have to look to the x-face that is one zone leftward of the left face of the zone of interest to obtain D x2− 0 . Figure 4 depicts this situation and shows the stencils. Using these two additional facial values, we can define two third order accurate stencils. The left-biased S 1 stencil consists of D x2− 0 , D x− 0 , D x+ 0 and the right-biased S 2 stencil consists of D x− 0 , D x+ 0 , D x2+ 0 . The stencils are also shown in Fig. 4. We can then write two quadratic polynomials, a left-biased polynomial d 1 (x) and a rightbiased polynomial d 2 (x) , as follows: Unlike finite-volume WENO reconstruction, the coefficients of these polynomials have to be tailored to match the facial values. Therefore, to specify d 1 (x) , we require Notice that our construction is such that both these polynomials will match D x− 0 and D x+ 0 at the left and right boundaries, respectively, of the zone being considered. Any convex combination of the two polynomials from the above two stencils will also match these facial values. Consequently, we can nonlinearly hybridize between the two polynomials in a WENO-like fashion. The coefficients of the polynomial on stencil S 1 are, therefore, given by and the coefficients of the polynomial on stencil S 2 are given by Because the polynomials in Eq. (14) use the same Legendre polynomial basis that were used in Balsara et al. [6], the smoothness indicators (which still consist of integrating the variation of the polynomial over the zone of interest) have identical expressions. Just for the sake of completeness, we provide them here If one wishes, one can arithmetically average d 1 (x) and d 2 (x) to obtain a third, centered stencil which can be given a higher weight just to prevent rapid stencil switching at this order. Note that at third order we have to accept three more polynomial terms in the full polynomial expansion to exactly match all the moments in the two opposing faces of the zone of interest. This completes our description of the third order WENO-ADP reconstruction.

r = 4 WENO-ADP
Say that the left and right x-faces of a coarse zone have up to piecewise-cubic transverse moments in the y-and z-directions so that we can write them as We now realize that to have all the volumetric moments up to fourth order, we only have to do moment-by-moment linear reconstruction in the x-direction for the quadratic and cubic terms in Eq. (18). For the linear terms, we have to use the r = 3 WENO-ADP reconstruction from the previous section, because we wish to obtain not just the xy-and xz-moments but also the x 2 y-and x 2 z-moments. We also need to reconstruct the constant, x-dependent, x 2 -dependent, and x 3 -dependent moments in the x-direction. Using Fig. 4, we now realize that we have three stencils S 1 , S 2 , and S 3 that rely on , respectively. The stencils are also shown in Fig. 4. Stencil S 1 is left-biased, S 2 is a central stencil, S 3 is a right-biased stencil, and all these stencils retain fourth order of accuracy. For i = 1, 2, 3 we can write the ith cubic polynomial corresponding to stencil S i as As before, we match the values at the faces of Fig. 4. For the polynomial on stencil S 1 we get the coefficients For the polynomial on stencil S 2 we get the coefficients For the polynomial on stencil S 3 we get the coefficients The ith smoothness indicator is given by Note that at fourth order we have to accept four more polynomial terms in the full polynomial expansion to exactly match all the moments in the two opposing faces of the zone of interest. Notice that five more terms is still substantially fewer than all the extra terms in Eq. (8), where we indeed match the facial moments as well as all the constraints. This completes our description of the fourth order WENO-ADP reconstruction.

r = 5 WENO-ADP
Say that the left and right x-faces of a coarse zone have up to piecewise-quartic transverse moments in the y-and z-directions so that we can write them as We now realize that to have all the volumetric moments up to fifth order, we only have to do moment-by-moment linear reconstruction in the x-direction for the cubic and quartic terms in Eq. (24). For the quadratic terms, we use the r = 3 WENO-ADP reconstruction, and for the linear terms we use the r = 4 WENO-ADP reconstruction. For the constant terms, we use the r = 5 WENO-ADP that is described for the rest of this subsection. Using Fig. 4, we now realize that we have four stencils S 1 , S 2 , S 3 , and S 4 that rely on , respectively. The stencils are also shown in Fig. 4. Stencil S 1 is extremely left-biased, stencil S 2 is partially leftbiased, stencil S 3 is partially right-biased, and stencil S 4 is extremely right-biased. For i = 1, 2, 3, 4 we can write the ith cubic polynomial corresponding to stencil S i as As before, we match the values at the faces of Fig. 4. For the polynomial on stencil S 1 we get the coefficients For the polynomial on stencil S 2 we get the coefficients For the polynomial on stencil S 3 we get the coefficients For the polynomial on stencil S 4 we get the coefficients The ith smoothness indicator is given by Note that at fifth order we have to accept five more polynomial terms in the full polynomial expansion to exactly match all the moments in the two opposing faces of the zone of interest. This completes our description of the fifth order WENO-ADP reconstruction. It also completed the section on WENO-ADP.

WENO-ADPT, a Touch-Up Procedure to Restore the Discrete Divergence on Fine Meshes with a Refinement Ratio of Two
The WENO-ADP reconstruction algorithm from the previous section can give us a high order representation of the vector field within a coarse zone. We first define a flux as an integration of the normal component of a vector field across an area. As a result, if a coarse zone is refined, the refinement procedure will generate new internal faces inside a coarse zone. To retain consistency on the coarse mesh, we do not touch the refined faces of a fine mesh that overlie a coarse mesh face. We only feel free to mildly touch-up the internal faces. In two-dimensions, Fig. 5 provides an example. Because the two-dimensional case is much easier to understand on our first encounter with the touch-up procedure, we present that first (by way of motivation) in Subsect. 4.1. In Subsect. 4.2, we present the three-dimensional case. (We are of course quick to add that the two-dimensional and threedimensional cases differ substantially in their complexity, and only the three-dimensional case is truly valuable to applications scientists.) All through this section, we restrict attention to refinement ratios of two. In the next section, we will address the case, where the refinement ratio can exceed two.

Touch-Up Procedure to Restore the Discrete Divergence on Two-Dimensional Meshes
This subsection is just meant to be an easy introduction to the topic with math that remains very simple because of the use of a two-dimensional example with a refinement ratio of two. Say that the coarse zone, shown with a solid boundary in Fig. 5, is refined with a refinement ratio of two. The newly defined fine mesh faces are shown with dashed lines.
The coarse mesh will have its top and bottom x-faces split into two with the result that we have the fluxes x t1 and x t2 at the top x-faces and the fluxes x b1 and x b2 at the bottom x-faces. Likewise, we have the fluxes  Designing to help us understand the touch-up procedure, WENO-ADPT, in its simplest two-dimensional form. It shows the subdivision of a coarse zone, shown by the solid lines, into four fine zones, whose boundaries are shown by the dashed lines. The bottom and top x-and y-faces of the coarse mesh have been subdivided into two. The area-integrated fluxes obtained from the vector field components at the bottom and top faces have subscripts "b" and "t" and they cannot be changed. The WENO-ADP reconstruction procedure for the vector field can be used in the interior of the zone to obtain the fluxes at the newly refined faces, and they are denoted with subscript "FV" and shown in red. These do not satisfy the divergence conditions, but they are higher order accurate. The volume-integrated charges that have to be matched in the four refined zones are shown by "Q 1 " to "Q 4 ". The "corrected" area-integrated fluxes, whose discrete divergence on the fine mesh exactly matches these charges, are shown in blue and have a subscript "c". The WENO-ADPT scheme is a procedure for obtaining the corrected fields in a fashion that keeps them as close as possible (in a constrained least squares sense) to the vector fields that have been obtained from the finite volume reconstruction refined faces in Fig. 5 are shown with dashed lines, and the WENO-ADP reconstruction can be used to endow them with the fluxes x , and y FV2 . These fluxes can have an order of accuracy that is given by the accuracy of the WENO-ADP scheme that we have chosen. We also have the finite volume, WENO-based, reconstructed charge density in the coarse zone. This can be volumetrically integrated to obtain the total charges Q 1 , Q 2 , Q 3 , and Q 4 in the four fine zones shown in Fig. 5. Note though that the fine mesh fluxes , and y FV2 , despite being high order accurate, will not satisfy an integrated divergence-preserving condition in each of the four fine zones in Fig. 5.
We will, however, have the following integrated divergence condition which holds over the coarse zone of Fig. 5: This shows us that of the four integrated divergence conditions that we can write for the four fine zones in Fig. 5, only three of them will be mutually independent. Now, let us define the fine mesh fluxes , and y c2 , which do indeed satisfy integrated divergence-preserving condition in each of the four fine zones in Fig , and y c2 . To fully solve the system, we supply a fourth parametric equation Equations (32) and (33) can be parametrically solved in terms of the parameter " ". Thus we get a parametrized solution to all the integrated divergence-preserving constraints given by Till the parameter " " is specified, the above equations are not very useful. They only illustrate the amount of freedom we have in the imposition of the divergence-preserving constraints.
It is most important to realize that we can indeed fully specify the parameter " " by minimizing the quadratic difference between the four divergence constraint-preserving , and y c2 and the four suitably high order fluxes x , and y FV2 . In other words, we wish to use the parameter " " to minimize Any reasonably good computer algebra system can easily give us the optimal parameter " " as With this optimal parameter " " inserted in Eq. (34), we get an analytically evaluated touch-up procedure on two-dimensional meshes with refinement ratios of two.
Since the flux was initially obtained by areal integration of the normal component of the vector field, the face-averaged normal components of the vector field that we actually need can be obtained by dividing the fluxes in Fig. 5 by their facial areas. Note too that because the divergence preserving-constraint is written in integral form, all the geometric dependences have been factored into the fluxes. As a result, we do not need to make a uniform subdivision of the coarse zone in Fig. 5. Furthermore, the boundaries of Fig. 5 do not need to be straight; they can indeed be curved and the mathematics remains unchanged as long as the fluxes are accurately computed. We see, therefore, that the touch-up procedure in our WENO-ADPT algorithm is indeed very general. This generality is retained in three dimensions.

Touch-Up Procedure to Restore the Discrete Divergence on Three-Dimensional Meshes
Now we pay attention to the touch-up procedure in the three-dimensional case. Consider Fig. 6a which shows a coarse zone and the nine x-faces that result when it is refined with a refinement ratio of two. Figures 6b and 6c show the same for the y-and z-faces. Figure 6d shows the labeling of the eight volume-integrated charges, Q 1 to Q 8 , in each of the eight fine zones. The sequence in which these fluxes and charges are labeled is very important especially for those seeking to make an implementation of the WENO-ADPT algorithm. Focusing on Fig. 6a, the four x-fluxes in the four refined top x-faces are given by , and x t4 . The four x-fluxes in the four refined bottom x-faces are given by , and x b4 . These eight fluxes remain unchanged in the touch-up procedure, because we want them to remain continuous across coarse mesh faces. The newly refined x-face is also shown and the WENO-ADP algorithm can be used to assign four x-fluxes to the refined faces and they are given by , and x FV4 . These four fluxes, along with their eight analogues from Figs. 6b and 6c, will not satisfy the integrated discrete divergence conditions in each of the eight refined zones. However, in Fig. 6a, we define four more , and x c4 . These four fluxes, along with their eight further analogues from Figs. 6b and 6c, will be made to satisfy the eight integrated discrete divergence conditions in the eight fine zones with the use of some free parameters. Therefore, the task at hand is to find the twelve constraint-satisfying fluxes, x , and z c4 which are closest to the twelve high order fluxes x , and z FV4 while simultaneously satisfying the integrated discrete divergence conditions in the eight fine zones. As in the previous subsection, this will be accomplished using the free parameters to minimize the squared difference between the two sets of fluxes.
Let us first establish the true degrees of freedom in the problem. Integrating over the coarse mesh, the divergence constraint becomes

3
This shows us that only seven of the eight divergence-preserving constraints in seven of the eight refined zones are indeed linearly independent. Therefore, on the fine mesh, we should account for seven independent divergence preserving constraints which we write in integrated form as Labeling of the y-fluxes at the bottom, "b", and top, "t", y-faces of a coarse zone. The newly refined y-face is also shown with the y-fluxes that are obtained from the WENO-ADP reconstruction shown in red and with subscripts "FV". The final y-fluxes, which have been corrected with the touch-up procedure, are shown in blue and with subscripts "c". c Labeling of the z-fluxes at the bottom, "b", and top, "t", z-faces of a coarse zone. The newly refined z-face is also shown with the z-fluxes that are obtained from the WENO-ADP reconstruction shown in red and with subscripts "FV". The final z-fluxes, which have been corrected with the touch-up procedure, are shown in blue and with subscripts "c". d Labeling of the eight volume-integrated charges in the eight fine zones that result from the refinement of the coarse zone Equation (38) gives us seven equations in twelve unknowns. To fully solve the system, we supply five more parametric equations as follows: Equations (38) and (39) can be parametrically solved in terms of the parameters " , , , , ". Thus we get a parametrized solution to all the integrated divergence-preserving constraints given by Till the parameters " , , , , " are specified, the above equations are not very useful. They only illustrate the amount of freedom we have in the imposition of the divergencepreserving constraints.
It is most important to realize that we can indeed fully specify the parameters " , , , , " by minimizing the quadratic difference between the twelve divergence constraintpreserving fluxes and the twelve suitably high order fluxes. In other words, we wish to use the parameters " , , , , " to minimize Any reasonably good computer algebra system can easily give us the optimal parameters " , , , , " as With these optimal parameters " , , , , " inserted in Eq. (30), we get an analytically evaluated touch-up procedure on three-dimensional meshes with refinement ratios of two. As before, the face-averaged normal components of the vector field that we actually need can be obtained by dividing the fluxes in Fig. 6 by their facial areas. Similarly, we do not need to make a uniform subdivision of the coarse zone in Fig. 6, the refined zones can be non-uniformly subdivided. In addition, as before, the boundaries of Fig. 6 do not need to be straight; they can indeed be curved and the mathematics remains unchanged as long as the fluxes are accurately computed. This concludes our description of the WENO-ADPT algorithm for prolonging divergence constraint-preserving vector fields in three dimensions when the refinement ratio is two.

WENO-ADPT, a Touch-Up Procedure to Restore the Discrete Divergence on Fine Meshes with Refinement Ratios Greater Than Two
Let us begin the discussion in this section by asking why this touch-up method, which is based on least squares minimization, works? For example, in the case, where the refinement ratio is two, we saw in the previous section that each refined zone will have 12 new facial variables in three dimensions. However, we also had seven mutually independent divergence conditions in the fine zones. As a result, we had five free parameters with which to bring the constraint-preserving fluxes in line with the fluxes that were obtained from the higher order WENO-ADP algorithm. Since 5/12 = 0.416 7, we see that there is considerable amount of freedom in the number of free parameters. As a result, the free parameters could always bring the values of the constraint-preserving fluxes very close to the values of the higher order fluxes. This is also illustrated in Fig. 7a which shows the number of constraints within each layer of the refined zones associated with refining a coarse zone by a refinement ratio of two. Figures 7b and 7c also show the number of constraints within each layer of the refined zones when refinement ratios of 3 and 4 are used. Figure 7 makes it easy to identify the number of constraints as a function of refinement ratio. Table 1 catalogues the number of newly introduced fine zone faces in the interior of a coarse zone as a function of increasing refinement ratio. It also catalogues the number of divergence constraints and the number of free parameters as a function of increasing refinement ratio. For all the refinement ratios tabulated in Table 1 we see that the ratio of the number of free parameters to the number of newly introduced facial variables in the interior of a newly refined coarse zone remains large. This means that there will always be sufficient number of free parameters with which to bring the constraint-preserving fluxes in line with the fluxes that were obtained from the higher order WENO-ADP algorithm. This is a very useful realization and enables us to understand why the touch-up algorithm (which is based on least squares minimization) works. Now that we have demonstrated the viability of the touch-up method with increasing refinement ratios, let us describe how the method is implemented. The newly refined fine mesh faces can be given an ordinal numbering, as shown in Fig. 6, and that numbering can be used for matrix assembly. Realize that the constraints are linear. From Fig. 7, we see that each refined zone contributes one linear constraint. The only exception is the last refined zone. As a result, assembling a matrix of linear constraints can be done in code, as can be the assembly of the corresponding right-hand side. Furthermore, the least squares problem is trivial. It just requires that the constraint-satisfying flux in each newly refined face (i.e., the flux with the subscript "c" in Fig. 6) should come as close as possible to the flux that is obtained from the suitably higher order WENO-ADP algorithm (i.e., the flux in the same face with subscript "FV" in Fig. 6). As a result, it is very easy to assemble the linear system for the least squares problem. Solving a constrained least squares (CLSQ) problem, therefore, emerges as the correct linear algebra-based approach for the touch-up procedure. The CLSQ method requires the assembly of a Karusch-Kuhn-Tucker (KKT) matrix and its corresponding right-hand side. For further details on the CLSQ method and the KKT matrix assembly, please see (Karusch [28], Kuhn and Tucker [29], Boyd and Vandenberghe [14]). If a fixed refinement ratio is used in all zones, and the refined mesh has zones of uniform size, the KKT matrix and its inversion have only to be done once at each refinement level. In each zone, we simply change the right-hand side and carry out a local, oncore, small matrix-vector multiplication to find the divergence constraint-satisfying fluxes. From these fluxes, the facial components of the desired constraint-preserving vector field are easily found, as described in the previous section. The algorithm is very efficient and embarrassingly parallel. This completes our description of the WENO-ADPT algorithm for prolonging divergence constraint-preserving vector fields in three dimensions on meshes with refinement ratios greater than two.

Implementation-Related Details
Here we provide stepwise implementation-related details for both the algorithms described in this work. In Subsect. 6.1, we describe the Algorithm from Sect. 2 at fourth order. Just to be specific, and also to show the higher order extension, in Subsect. 6.2 we describe the Algorithms from Sects. 3 to 5 at fifth order. We also discuss ways to make the implementation more efficient.

Implementation of the Algorithm from Sect. 2
(i) Using the results from Supplement A.2, reconstruct the facial modes for those coarse zone faces that do not abut a fine mesh. For those faces that do abut a fine mesh, use Eq. (7) to obtain the facial modes. This should fully specify all the facial modes in Eqs.
(2) to (4). (ii) Using Eq. (5b), obtain the charge density within each coarse zone. (iii) Using the results from Supplements A.2 and A.4, obtain all the higher moments of the charge density. (iv) At this point, all the facial modes and the volumetric charge density in the coarse zone have been specified up to fourth order of accuracy. Therefore, using the results from Supplement B we can fully specify all the constraint-preserving modes of the vector field in Eqs. (8) to (10). (v) Once this is specified, the field components on the fine mesh can be obtained by areal integration and area-weighted averaging. By design, this assignment will be fourth order accurate and globally constraint-preserving.

Implementation of the WENO-ADPT Algorithm from Sects. 3, 4, and 5
(i) Using the results from Supplement A.3, reconstruct the facial modes for those coarse zone faces that do not abut a fine mesh. For those faces that do abut a fine mesh, use Eq. (7) (or any suitable higher order extension) to obtain the facial modes. This Use these values to obtain fluxes. This will be a fifth order accurate prolongation. These are the fluxes with subscripts "b" and "t" in Fig. 6. (vii) Now using the reconstructed WENO-ADP modes, which indeed extend into the volumes of the coarse mesh, we prolong the components of the vector field to those faces of the fine mesh that do not coincide with the coarse mesh faces. This will be a fifth order accurate prolongation, but at this stage of the game it will not be a globally divergence constraint-preserving prolongation. Use these values to obtain fluxes. These are the fluxes with subscripts "FV" in Fig. 6. (viii) For a refinement ratio of two, use Eqs. (40), (41), and (42) to obtain the fifth order accurate, divergence constraint-preserving fluxes. These are the fluxes with subscripts "c" in Fig. 6. From the fluxes, retrieve the vector field components within the faces of the refined mesh. These vector field components will indeed be fifth order accurate and divergence constraint-preserving. (ix) For a refinement ratio greater than two, use the CLSQ process from Sect. 5 to do analogously to the previous step. From the fluxes, retrieve the vector field components within the faces of the refined mesh. These vector field components will indeed be fifth order accurate and divergence constraint-preserving. 3 × (r 3r 2 ) r 3 -1 2 × r 3 -3 × r 2 + 1 In the limit of r → ∞, the ratio → 2/3 = 0.667

Opportunities for Even More Efficient WENO Implementation
In Jiang and Shu [27] and Balsara and Shu [9] a class of WENO algorithms were presented and shown to extend to all orders in finite difference context. In Balsara et al. [6] we showed that even more compact expressions could be derived for the stencils along with very compact, perfect-squares expressions for the smoothness indicators. The supplementary materials in this paper show that even in the context of finite volume WENO, compact expressions can be derived for the stencils along with very compact, perfect-squares expressions for the smoothness indicators. These innovations contribute to the efficient implementation of WENO schemes. In this paper, we have documented these innovations for all possible stencils so that the users of WENO schemes have all the options available to them. However, various practitioners have suggested even more efficient ways of assembling the stencils so that even fewer zones are used for the reconstruction. It may be valuable to very briefly document some of the opportunities for using fewer zones in stencil-construction here.
In Balsara and Shu [9] itself it was realized that the higher order WENO variants could use some extra monotonicity preservation. The central WENO (CWENO), WENO-ZQ, and WENO with adaptive order (WENO-AO) schemes (Cravero and Semplice [20], Cravero et al. [19], Zhu and Qiu [42], Balsara et al. [6]) also derived their robustness from the inclusion of smaller stencils in addition to the use of one large central stencil that preserved the order property. CWENO and WENO-ZQ used piecewise linear reconstruction as their choice of smallest stencil. Using third-order WENO as the smaller stencil, Balsara et al. [6] were able to retain the prospect that local extrema could be somewhat represented on the mesh. See also Zhu and Shu [43] as another alternative. (Here we have focused exclusively on structured meshes, because that is the focus of this paper.) All these methods make the overall stencil operations much simpler, and the overall stencil much smaller, contributing to the efficiency of WENO implementation.
All the methods in the above paragraph are agnostic to the underlying PDE. In other words, regardless of the solution characteristics of the PDE system, they try to achieve their efficiencies by making the stencil operations more efficient. If one can factor in some knowledge of the underlying PDE, then one can achieve even greater efficiencies (Li and Ren [30]). In such approaches, one can dispense with the WENO stencil operations altogether and use just the central stencil in parts of the domain, where the solution is known to be indifferent to the need for nonlinear hybridization. Only a few troubled zones would need the WENO nonlinear hybridization. However, such methods require practitioners with a superb and intuitive feel for the underlying PDE.

Description of the Update Algorithm on a Composite Mesh
In Balsara [1], a second-order divergence-free AMR algorithm was presented which used a one-step, Heun-type, predictor-corrector scheme for its update. The availability of ADER methods makes it possible to use predictor steps that have better than second-order accuracy in time (Dumbser et al. [23]), and the availability of the reconstruction methods described here makes it possible to have higher order accuracy for divergence-preserving PDEs in space. As a result, working in the context of a single-stage higher order predictor-corrector scheme, it is possible to upgrade the roadmap from Balsara and Shu [9] so as to make it suitable for any order of accuracy. We assume that the underlying PDE has one or more divergence-constraints. Here we describe a divergence-constraint preserving algorithm for a two-level composite mesh. However, the reader will easily see that the method recursively extends to any composite/adaptive mesh hierarchy. (For our purposes, a composite mesh just refers to the fact that we view the time update on the fine and coarse mesh as a single entity.) We assume a refinement ratio of two just to keep the narrative simple, but once the ideas are understood they naturally extend to other refinement ratios. To keep everything very simple, we assume a Faraday's law update equation ( t + ∇ × = 0 ), where the facially collocated magnetic induction (B) is updated by the circulation of the edge-collocated electric field (E), Figs. 1b and 1c.
AMR meshes have to be processed with optimal spatial and temporal efficiency. On a two-level mesh, the coarse mesh goes from time t n to time t n+1 = t n + Δt by taking a single timestep of size Δt in time. Because of Courant number restrictions, the fine mesh can only go from time t n to time t n+1∕2 = t n + Δt∕2 during its first timestep and from time t n+1∕2 to time t n+1 = t n+1∕2 + Δt∕2 during its second timestep. As a result, the fine and coarse meshes time-synchronize at times t n and t n+1 . We wish to describe an algorithm that is so simple that much of the computational complexity is expended on updating single (coarse or fine) meshes with suitable amounts of ghost zones. To that end, focus on Fig. 8. It shows the simplest schematic diagram of a two-level adaptive mesh. The solid blue lines show a coarse mesh. The ghost zones for that mesh are not shown. The solid red lines show a fine mesh. The ghost zones for that fine mesh are indeed shown with dashed red lines in Fig. 8. We show only two layers of ghost zones, suitable for a second-order scheme, but the number of layers of ghost zones can be increased with increasing order.
In AMR, it is always assumed that one has "proper nesting" between any fine mesh and its parent coarse mesh. In other words, it is assumed that the fine mesh is either fully nested within a coarse mesh, as shown in Fig. 8, or that the fine mesh coincides with a physical boundary (in which case boundary conditions apply to the part of the fine mesh that coincides with the boundary). If a fine mesh is too close to a boundary, but does not coincide with a boundary then it is a good idea to make it properly nested by extending it all the way to the boundary! At time t n we can make a divergence-preserving reconstruction on the coarse mesh. If the fine mesh has been newly built, the coarse mesh (after reconstruction) has all the information for initializing the fine mesh and its ghost zones in divergencepreserving fashion. Even if the fine mesh is a pre-existing patch of fine mesh, its ghost zones should be filled in via prolongation of the reconstructed divergence-preserving solution from the coarse mesh at time t n . As a result, the fine mesh has everything it needs to take the first of its two timesteps.
The coarse mesh is always updated in time before the fine mesh. Therefore, let us assume that a timestep has been taken on the coarse mesh, bringing its solution to a time of t n+1 = t n + Δt . In the course of taking that timestep, an ADER predictor step is effected on the coarse mesh. The ADER formulations that we prefer are given in Dumbser et al. [23] or the more recent ADER formulation in Balsara et al. [11]. Either formulation gives us a full and highly accurate space-time representation of the solution of the PDE within each coarse mesh zone. Now we make the fine mesh take its first timestep from time t n to time t n+1∕2 = t n + Δt∕2 . At time t n+1∕2 , realize that the ghost zones of the fine mesh in Fig. 8 will need to be updated in time-accurate fashion. This is where the space-time accurate ADER predictor step on the coarse mesh comes in. The ADER space-time reconstruction is stored on the coarse mesh at fine-coarse interfaces and it allows us to provide good values to the fine mesh ghost zones at time t n+1∕2 . (Please understand that the ADER solution started from a divergence-preserving solution so it will be very accurate and quite close to divergence-preserving. Therefore, it can be used to fill in ghost zones, and only the ghost zones, at time t n+1∕2 but it can never be used to re-initialize a fine mesh.) The fine mesh now has everything it needs to take its second timestep from time t n+1∕2 to time t n+1 = t n+1∕2 + Δt∕2 . A schematic diagram of the narrative in this paragraph is shown in Fig. 9. The thick black horizontal lines indicate the time levels in the update of the composite mesh. The yellow ovals in Fig. 9 show the prolongation steps, where the fine mesh ghost zones are filled with information that is available from the coarse mesh. The upwardpointing blue arrow indicates one coarse mesh timestep. The two upward-pointing red arrows indicate that we are taking two fine mesh timesteps for one coarse mesh timestep.
In principle, it would seem at this point in the narrative that the fine and coarse meshes are time-synchronized at time t n+1 , but this synchronization is far from perfect and we have to consider a further point of detail. Consider Fig. 10 of this paper (which is adapted from Balsara [1]) which shows an exploded view of a fine-coarse interface in the mesh hierarchy. Our philosophy in AMR is that the finest mesh has the most accurate solution and its solution should be represented as much a possible on the composite mesh. Therefore, at the fine-coarse interface, we have a dilemma, because the capital B-fields of the coarse mesh have been updated with the capital E-fields of the coarse mesh, see Fig. 10. However, the lower-case b-fields of the fine mesh have been updated with the lower-case e-fields of the fine mesh, see Fig. 10. We see from the exploded view in Fig. 10 that the upper and lower case fields coincide at the fine-coarse interface. Using our previously stated philosophy, we should give primacy to the lower case, i.e., fine mesh, variables. The overlapping facially collocated and edge-collocated variables at the fine-coarse interface shown in Fig. 10 can be restored by a process of defect correction. Defect-correction in the AMR context means that we have to undertake a special sequence of operations at time t n+1 to restore consistency at fine-coarse interfaces. Because of our philosophy, it is the coarse mesh variables that will have to be given corrections to make the coarse mesh solution consistent with our (more accurate) fine mesh variables. Defect correction is also schematically shown by the purple triangles and the brown rectangle in Fig. 9. We explain some more details about defect correction in the next two paragraphs; but the best description is still the one in Sect. 5 of Balsara [1].
Say that the coarse mesh has zones of size Δx , Δy , and Δz in the x-, y-, and z-directions and say that the fine mesh has zones that are half as small. On the coarse mesh, say that the one-step predictor-corrector formulation produces electric fields with a superscript of "n + 1/2". We can write the update from time t n to time t n+1 on the coarse mesh as Please see Fig. 10 , and E n+1∕2 2,y are higher order time-averaged quantities and can be better than second-order accurate. We take just one instance of the two update steps on the fine mesh. Let the electric fields that take the fine mesh solution from time t n to time t n+1∕2 be denoted with a superscript "n + 1/4" and let the electric fields that take the fine mesh solution from time t n+1∕2 to time t n+1 be denoted by a superscript "n + 3/4". The two timestep equations for the update of b n 1,x can be written as (Equations that are analogous to the above two equations can be written for b Fig. 10.) As before, the electric fields used in the updates in Eqs. (44) and (45) are higher order time-averaged quantities and can be better than secondorder accurate. Now notice that we start the timestep on the composite mesh with and we would like to end the timestep on the composite mesh with Note though that Eqs. (43) to (45) do not guarantee that Eq. (47) will hold at the end of a composite mesh update. This is because we cannot guarantee the following equalities for the electric fields on coarse and fine meshes at the locations, where the electric fields coincide in Fig. 10. The desired equalities are However, in practical computation, the above equalities are never exactly attained. This is where we have to go back to our philosophy from the previous paragraph-wherever   The ghost zones for that mesh are not shown. The solid red lines show a fine mesh. The ghost zones for that fine mesh are indeed shown with dashed red lines. When the fine mesh is initialized, the solution is prolonged from the blue coarse mesh to all locations of the red fine mesh, including the ghost zones. Whenever the two meshes synchronize in time, the prolongation is also used to initialize the ghost zones of the fine mesh they coincide, the fine mesh variables (because they are more accurate) should be given precedence over the coarse mesh variables! This is what the defect correction step is intended to accomplish. In another way of looking at it, Eq. (47) would indeed have been exact if the update of the coarse mesh magnetic induction had been Note that we do not have the electric fields on the right-hand side of Eq. (49) when we begin to initially take the timestep on the coarse mesh. As a result, out of practical necessity, we have to initially take the coarse mesh timestep using Eq. (43). However, please also observe that by the time the composite mesh is time-synchronized at time t n+1 , all the electric fields on the right-hand side of Eq. (49) have been evaluated. As a result, if we can intelligently store them in a suitable way, we can, at the moment of time-synchronization, effectively achieve Eq. (49). This is the idea behind defect correction. In the next paragraph we explain practical details of defect correction.
As a matter of practical necessity, it is easiest to store the defects in the electric field on the coarse mesh. Specifically, we store them at the same locations that correspond to the  fine-coarse interface on the coarse mesh. Now please see Fig. 9. We realize that after the coarse mesh has completed its timestep, we retain only the (negative of the) coarse mesh electric fields that correspond to the fine-coarse interface. As the fine mesh undergoes its two updates, the fine mesh electric fields that coincide with that interface are added, with a suitable coefficient of 1/4, to the (negative of those same) coarse mesh electric fields. As a result, the defect is assembled on the coarse mesh as the timestep proceeds on the composite mesh. When the fine and coarse meshes are time-synchronized again at time t n+1 we can apply the defect correction to the coarse mesh magnetic induction that has already been updated using Eq. . Indeed, the defect correction gives primacy to the fine mesh electric fields over the coarse mesh electric fields at  in Fig. 10. The curly bracketed terms in Eqs. (50), (51), and (52) are indeed the defects that are assembled over the course of the timestep on the composite mesh. Please see the purple triangles and the brown rectangle in Fig. 9. The purple triangles indicate the act of assembling the defects as the timestep proceeds on the composite mesh. In other words, the purple triangles show how the curly brackets in the above three equations are assembled and stored at the fine-coarse interface on the coarse mesh. The brown rectangle indicates the act of enforcing the defect correction on the coarse mesh as shown in Eqs. (50), (51), and (52). We hope the reader has realized by now that the process of defect correction simply consists of assembling the defects in the electric field at the fine-coarse interfaces, and nowhere else. We then apply one time-update step on the coarse mesh, where the defects in the electric field are used instead of the actual electric fields. In other words, if the user has written one nice subroutine for the update step, that same subroutine can be reused for the defect-correction step, where we use the defects in the electric fields rather than the electric fields themselves. Recall too that the defect is zero at all locations except at the fine-coarse interfaces.
The defect correction, described in the above two paragraphs, restores consistency at the fine-coarse interface. It also makes the entire composite mesh divergence-preserving (or divergence-free depending on what is desired). Any inability to retain this consistency will show up as instabilities in the code or as a localized build-up of divergence at the finecoarse interface. Because these deficiencies are small over a single composite timestep, but pernicious over time, they will show up as late-time instabilities in any code that is bereft of defect-correction that processes a composite mesh with timestep sub-cycling. Now please recall once more our philosophy of giving primacy to the fine mesh variables. It necessarily means that when the fine and coarse meshes synchronize at time t n+1 , all the fine mesh variables should be given primacy on the composite mesh. This is done via the restriction step. It is shown by the green pentagon in Fig. 9. In this step, we simply take all fine mesh data and average it suitably to coinciding locations on the coarse mesh. To take an example from Fig. 8, when the fine and coarse meshes are time-synchronized, the red faces on the fine mesh should be used to reset the coinciding blue faces on the coarse mesh. This reset will look exactly like Eq. (47). This completes our description of the restriction step. At this point we have described all the steps in Fig. 9, so our description of composite mesh timestep-processing is complete.
As a last, precautionary note, it is worth pointing out one more important point. There will be times when one wants to extract the composite mesh information on a global fine mesh so that one can analyze the results further. This can happen, for instance, when one is doing accuracy analysis. On such occasions, please be mindful of one more point of detail at fine-coarse interfaces. At those interfaces, the facial reconstruction on the coarse mesh should be made consistent with the facial reconstruction on the overlapping faces of the fine mesh. Equation (7) provides an example of how this is done.

Results from Prolongation
We can see that the results described in Sect. 2 are analytically exact at fourth order. Furthermore, the mathematical construction works for divergence-free and divergence-preserving vector fields with refinement ratios of two. Therefore, in Subsect. 8.1, we demonstrate that divergence-free prolongation from a coarse mesh to a fine mesh with refinement ratio two that uses Sect. 2 indeed works. In Subsect. 8.2 we demonstrate that divergencepreserving prolongation from a coarse mesh to a fine mesh with refinement ratio two that uses Sect. 2 indeed works. Subsections 8.1 and 8.2 show that a simple approach is available for the most basic need, which is to provide a higher order accurate divergence-free and divergence-preserving prolongation strategy that can work with refinement ratio of two. However, the results in those sections also provide us with a point of reference for the subsections that are to follow which show the versatility of the WENO-ADPT algorithm working at all orders and with different refinement ratios. We display refinement ratio of two, which exercises the results from Sect. 4, and refinement ratio of three, which exercises the results from Sect. 5. Subsections 8.3 and 8.4 show results from the WENO-ADPT algorithm for carrying out divergence-free and divergence-preserving prolongation, respectively, when refinement ratios of two are used. Subsections 8.5 and 8.6 show results from the WENO-ADPT algorithm for carrying out divergence-free and divergence-preserving prolongation, respectively, when refinement ratios of three are used. Subsection 8.7 shows the utility of using the finite volume WENO developed here for the prolongation of a scalar variable, or in fact, any vector field that does not have an associated constraint.

Divergence-Free Fourth-Order Prolongation with Refinement Ratio Two, Using
Sect. 2 The best way to generate a divergence-free field is to start with a vector potential and take its curl, ∇ × . The resulting vector field will be divergence-free. If we are starting with an analytically-specified vector potential, we can even ask a computer algebra system to provide the facially averaged normal components of ∇ × , so that they can be assigned to the faces of the coarse mesh that we start with. Here we start with the vector potential: The coarse mesh has an extent of [−0.5, 0.5] 3 . A fine mesh is then specified with a refinement ratio of two. The fine mesh zones uniformly sub-divide the coarse mesh so that each coarse zone is divided into eight equally sized fine mesh zones. The fine mesh always has eight times as many zones as the coarse mesh. Because the vector field can be analytically specified on the fine mesh too, we can use that information to extract the L 1 and L ∞ (53) ⎧ ⎪ ⎨ ⎪ ⎩ A x (x, y, z) = sin(xπ) cos(yπ) cos(3zπ); A y (x, y, z) = cos(yπ) cos(4xπ) cos(zπ); A z (x, y, z) = sin(zπ) cos(xπ) cos(2yπ) .
norms of the errors on the fine mesh. We can also use the fine mesh to verify that the prolongation has been divergence-free. Table 2 shows the results. We see that the x-component of the vector field has been prolonged with fourth order of accuracy. The other components also have the same order of accuracy, but to save space, they are not shown in the table.
Only the accuracy of the fine faces that do not coincide with the coarse faces is shown, because we have found that the coinciding fine mesh faces have superlative accuracies. The maximum absolute value of the divergence of the vector field on the fine mesh remains close to machine precision. The L ∞ norm in Table 2 only measures the maximum pointwise deviation and for the WENO schemes described here. At fourth order, the methods from Sect. 2 carry out a lot of float point evaluations. Therefore, the L ∞ norm is not as well-regulated as one would desire. However, the important point is that even with a very basic WENO the L 1 error, nevertheless, reaches fourth order. We will see in the ensuing discussions that the methods from Sects. 3 and 4 perform better in that regard.

Divergence-Preserving Fourth-Order Prolongation with Refinement Ratio Two, Using Sect. 2
When considering a vector field that is not divergence-free, one can start with any vector field. This is because any randomly chosen vector field will necessarily generate some amount of divergence. We choose the vector field given by The coarse mesh has an extent of [−0.5, 0.5] 3 . A fine mesh is then specified with a refinement ratio of two. The fine mesh zones uniformly sub-divide the coarse mesh so that each coarse zone is divided into eight equally sized fine mesh zones. The fine mesh always has eight times as many zones as the coarse mesh. Because the vector field can be analytically specified on the fine mesh too, we can use that information to extract the L 1 and L ∞ norms of the errors on the fine mesh. We can also use the fine mesh to verify that the prolongation has been divergence-preserving. The L 1 norm for the charge density is also shown, because we wish to show that the volumetric prolongation of the charge density using the finite volume WENO developed in Supplement A also satisfies an order property. Table 3 shows the results. We see that the x-component of the vector field has been (54) y, z) = 3x sin(xπ) cos(yπ) cos(4zπ); D y (x, y, z) = 2y cos(yπ) cos(2xπ) cos(zπ); D z (x, y, z) = 5z sin(zπ) cos(xπ) cos(3yπ) . prolonged with fourth order of accuracy. The other components also have the same order of accuracy, but to save space, they are not shown in the table. Only the accuracy of the fine faces that do not coincide with the coarse faces is shown, because we have found that the coinciding fine mesh faces have superlative accuracies. The charge density has also been prolonged with fourth order of accuracy. The maximum absolute value of the residual of the divergence constraint on the fine mesh remains close to machine precision. However, please note that as we go to more refined meshes, the divergence is evaluated via division by a progressively smaller zone size. As a result, there is a slow secular increase in |∇ ⋅ − | max in Table 3 (and tables like it) as the mesh resolution is increased. This is nothing to be alarmed of because it is purely a result of division by a progressively smaller number (the zone size). In Sect. 9 we will show that the undivided divergence remains nicely bounded with increasing refinement, and over very large numbers of timesteps.

Divergence-Free Prolongation at Second to Fifth Orders with Refinement Ratio Two, Using WENO-ADPT Algorithm from Sects. 3 and 4
Here again, we use the setup from Eq. (53). As before, each coarse zone is divided into eight equally sized fine mesh zones. We show the L 1 and L ∞ norms of the errors of the x-component of the vector field evaluated at the non-overlapping faces of the fine mesh. We also use the fine mesh to verify that the prolongation is initially not divergence-free, but we show how the touch-up procedure eventually makes it so. Table 4 shows the results when the second to fifth order WENO-ADP algorithms from Sect. 3 are applied, i.e., when there is no touchup. We see that the vector field preserves order of accuracy, illustrating that the WENO-ADP reconstruction algorithms do meet their design accuracies. However, the vector field is not divergence-free. We can see, nevertheless, that the divergence (evaluated in differential form) is one order less accurate than the design of the scheme used. This shows us that the error build-up in the divergence is well controlled. Table 5 shows the same results as Table 4, but after the application of the touch-up procedure from Sect. 4. We see that the order of accuracy for the vector field is intact. If anything, we observe that the error in the vector fields in Table 5 is somewhat lower than the error in the vector fields in Table 4. This shows that the imposition of the constraint actually reduces the error. We also see that the maximum absolute value of the divergence of the vector field on the fine mesh has been brought close to machine precision. This shows that the WENO-ADPT algorithm has worked at all the orders at which we tested it. We can also compare the fourth order WENO-ADPT results from Table 5 to the results  from Table 2. We find that the error and order of accuracy are overall competitive. However, we see that the error in Table 2 is somewhat lower than the error in Table 5, especially on very coarse meshes. This is comprehensible, because the WENO-ADPT uses multiple reconstruction steps, whereas the results in Table 2 were obtained with just one reconstruction step. Therefore, we see that reconstruction causes some degradation of error, especially on very coarse meshes. On finer meshes, the errors in the vector field from Tables 2 and 5 become comparable. Another interesting insight can be gained by looking at the column of the divergence. In this regard, we see that Table 5 has superior results compared to Table 2. This also makes sense, because the least squares minimization procedure that was used in Table 5 directly imposes the constraints on the fine mesh. Notice that the fourth order divergence-preserving reconstruction from Sect. 2 and Supplement B uses a lot of float point operations, and the rounding errors in those operations are reflected in Table 4 Results (before the touch-up procedure) of the prolongation of the x-component of the vector field at orders ranging from second to fifth. Here we show only the L 1 and L ∞ norms of the errors in the x-component of the vector field. The maximum of the absolute value of the divergence of the vector field on the fine mesh zones is also shown; along with the accuracy of the constraint. A refinement ratio of two was used the slightly degraded constraint-preservation in Table 2. In summary, the algorithm presented in Sect. 2 is slightly more accurate (especially on very coarse meshes) because of the smaller number of reconstruction steps, but as a counterpoise, the WENO-ADPT algorithm from Sects. 3 and 4 is better at preserving the constraints.

Divergence-Preserving Prolongation at Second to Fifth Orders with Refinement Ratio Two, Using WENO-ADPT Algorithm from Sects. 3 and 4
For this divergence-preserving demonstration, we again pick the setup from Eq. (54). The fine mesh is a uniform subdivision of the coarse mesh. We also use the fine mesh to verify that the prolongation is initially not divergence-preserving, but we show how the touch-up procedure eventually makes it so. Table 6 shows the results when the second to fifth order WENO-ADP algorithms from Sect. 3 are applied, i.e., when there is no touch-up. The algorithm is seen to preserve order of accuracy, illustrating that the WENO-ADP reconstruction algorithms do meet their design accuracies. However, the vector field is not divergence-preserving. Even so, we can observe that the maximal residual of the divergence Table 5 Results (after the touch-up procedure) of the prolongation of the x-component of the vector field at orders ranging from second to fifth. Here we show only the L 1 and L ∞ norms of the errors in the x-component of the vector field. The maximum of the absolute value of the divergence of the vector field on the fine mesh zones is also shown. A refinement ratio of two was used condition is only one order less accurate than the design of the scheme used. This shows us that the error build-up in the divergence is well controlled. Table 7 shows the same results as Table 6, but after the application of the touch-up procedure from Sect. 4. We see that the order of accuracy of the vector field is still preserved, while the divergence constraint is driven down to machine accuracy. This shows that the WENO-ADPT algorithm has worked at all the orders at which we tested it. We also see that the charge density retains its design accuracy. Since the prolongation of the charge density used the finite volume WENO algorithm from Supplement A, it shows that WENO algorithms are excellent vehicles for the prolongation of scalar fields or vector fields that do not have any particular constraint. We can also compare the fourth order WENO-ADPT results from Table 7 to the results from Table 3. We find that the error and order of accuracy are overall competitive. We again see that Table 3 shows more accurate results on coarser meshes than Table 7. However, the Table 6 Results (before the touch-up procedure) of the prolongation of the x-component of the vector field at orders ranging from second to fifth. Here we show only the L 1 and L ∞ norms of the errors in the x-component of the vector field. The maximum absolute value of the residual of the divergence constraint on the fine mesh is also shown; along with the accuracy of the constraint. A refinement ratio of two was used Coarse mesh resolution L 1 error for Bx L 1 accuracy for Bx L ∞ error for Bx L ∞ accuracy for Bx constraint-preservation is superior in Table 7. The explanations for this behavior that were given in the previous subsection carry over to this subsection.

Divergence-Free Prolongation at Second to Fifth Orders with Refinement Ratio Three, Using WENO-ADPT Algorithm from Sects. 3 and 5
This subsection mirrors Subsect. 7.3, with the only difference that a refinement ratio of three was used. Table 8 shows the results of applying the WENO-ADP algorithm. Table 9 shows the results of applying the touch-up procedure, i.e., the full WENO-ADPT algorithm. As before, even when the refinement ratio is larger than two, we see that the method works, which is a nice illustration of its generality.

Divergence-Preserving Prolongation at Second to Fifth Orders with Refinement Ratio Three, Using WENO-ADPT Algorithm from Sects. 3 and 5
This subsection mirrors Subsect. 7.4, with the only difference that a refinement ratio of three was used. Table 10 shows the results of applying the WENO-ADP algorithm. Table 11 shows the results of applying the touch-up procedure, i.e., the full WENO-ADPT algorithm. As before, even when the refinement ratio is larger than two, and even when the vector field is not divergence-free but indeed divergence-preserving, we see that the method works. Again, this is a nice illustration of its generality.

Accuracy of the Finite Volume WENO Scheme; When Used for Prolongation of Scalars
In addition to the prolongation of the constrained vector fields, we have also presented a very efficient finite volume WENO scheme. That scheme was used all through for the prolongation of the charge density. Therefore, we can also document the accuracy of the finite volume WENO method when it is used as an algorithm for prolonging scalars (and also unconstrained vector fields). Table 12 shows the results. We see that the finite volume WENO is a very effective algorithm for the prolongation of unconstrained variables on AMR hierarchies.

Time-Dependent AMR Calculation Showing a Resolution of the Late-Time Instability in CED
In Sect. 7, and via Figs. 8, 9, and 10, we have presented a set of algorithmic ingredients that permit temporally sub-cycled simulations on AMR hierarchies. When FDTD methods are used on composite meshes, the methods become susceptible to a late-time instability Table 9 Results (after the touch-up procedure) of the prolongation of the x-component of the vector field at orders ranging from second to fifth. Here we show only the L 1 and L ∞ norms of the errors in the x-component of the vector field. The maximum of the absolute value of the divergence of the vector field on the fine mesh zones is also shown. A refinement ratio of three was used (Liu and Sarris [32,33] and references therein). Such a late-time instability shows up in long-running CED calculations that use FDTD on composite meshes. The instability arises when multiple waves impinge on the interface between a coarse mesh and a fine mesh. The instability is exacerbated when the fine-coarse interface does not move. The instability reveals itself via unphysical growth of the solution at the fine-coarse interface. For that reason, in this section we construct a test problem with those characteristics and show that the methods developed here are free of the above-mentioned deficiencies. Our further goal is to construct a test problem, drawn from CED, that can be used to demonstrate the order property on composite meshes. As a result, in this section we present a time-dependent two-level CED test problem. We show that the solution using our methods is free of latetime instability and also preserves its designed second, third, and fourth orders of accuracy. The test problem was constructed with certain priorities. First, we wanted a simple composite mesh structure that is easy to implement and run in a serial setting. Second, we wanted it to clearly show the accuracy of the method on an entirely non-trivial problem for which there is no analytical solution. For this reason, we would have to run the problem on a very high resolution uniform mesh to obtain a reference solution which can then be used to extract a numerically motivated order of accuracy. Third, we wanted a problem which shows impinging of multiple waves from the coarse mesh onto the fine mesh so that we can demonstrate that the resulting method is free of late-time instability. Because FV-type schemes only show their asymptotic order of accuracy at high resolutions, we had to pick a two-dimensional problem.
The equations to be solved are Maxwell's equations for the propagation of waves in a dielectric medium. For Maxwell's equations, the primal variables are the electric displacement and the magnetic induction which evolve according to the curl of the magnetic field and the curl of the electric field . The constitutive relations between these vector Table 11 Results (after the touch-up procedure) of the prolongation of the x-component of the vector field at orders ranging from second to fifth. Here we show only the L 1 and L ∞ norms of the errors in the x-component of the vector field. The maximum absolute value of the residual of the divergence constraint on the fine mesh is also shown. A refinement ratio of three was used fields are = and = , where is a 3 × 3 permittivity tensor and is a 3 × 3 permeability tensor. For most materials these tensors tend to have a constant value along the diagonal and zero off-diagonal elements. Because of the curl-type evolution, and satisfy involution constraints that are a natural part of Maxwell's equations. Maxwell's equations can be written in MKS units as In this problem the current density is e = 0 and the charge density is e = 0 . The first equation in the above set is referred to as the generalized Ampere law; the second equation in the above set is referred to as Faraday's law; the third equation is called Gauss' law and the fourth equation embodies the fact that magnetic monopoles are absent. The entire problem is specified in MKS units in the xy-plane. Here we simplify the problem even further by taking = diag r 0 , r 0 , r 0 and = diag 0 , 0 , 0 . Here 0 = 8.85 and 0 = 4π × 10 −7 T ⋅ m ⋅ A −1 are the permittivity and permeability of free space. The dimensionless r is the relative permittivity and can vary in different media. The speed of light is given by c ≡ 1 � √ r 0 0 . Figure 11 shows the variation of the relative permittivity in the two-dimensional domain spanning (x, y) ∈ [−14, 14] × [ −7, 7] , where the distances are in meters. It shows a cylinder with enhanced permittivity, where the relative permittivity is specified along the cylinder's radius by The relative permittivity has been given a taper so that the same permittivity variation can be consistently assigned to meshes with different resolutions; this is essential for carrying out an accuracy analysis. When the mesh refinement is used, we use a refined patch that is centered on the cylinder and extends from [−7, 7] × [−3.5, 3.5] , again with square zones that have linear dimension in meters that is half that of the coarse mesh. The rectangle with solid lines inside Fig. 11 shows the extent of the refined patch. The cylinder is impacted from the left by a Gaussian pulse of electromagnetic radiation, where the pulse profile contains several wavelengths worth of electric and magnetic fields. The pulse is centered at (a, b) = (−10.5, 0) and the majority of the radiation energy in the pulse moves in the positive x-direction at t = 0. All the boundaries of the global mesh are continuative outflow. The waves are TM z waves, where the electric displacement is in the xy-plane and the magnetic induction is in the z-direction, i.e., perpendicular to the xy-plane. The pulse is initialized with a magnetic vector potential given by and an electric vector potential given by The resulting magnetic induction field ≡ ∇ × and electric displacement field ≡ c 0 (∇ × ) are discretized on the mesh using a discrete version of Stokes law. The simulation proceeds for (17 m)∕ 3 × 10 8 m/s ≈ 5.7 × 10 −8 s , by which time the incoming Gaussian pulse has scattered strongly off the dielectric cylinder. The problem is run with a CFL of 0.4. On finer meshes, this stopping time corresponds to many thousands of timesteps, thus enabling us to probe the late-time behavior of this problem. Figure 12 shows the solution variables Dx, Dy, and Bz for the test problem, as shown in Fig. 11. We used a second-order scheme with a 2 880 × 1 440 zone mesh. At that resolution, any further change in the solution becomes indistinguishable by eye. The solution on this mesh was used as a reference solution for all our second-order runs. Similarly high resolution third and fourth order simulations were used as reference simulations for our third and fourth order tests. at which we carry out accuracy analysis. The black circle shows the extent of the disk with a radius of unity. Table 13 shows results from the second-order accurate runs. We see the L 1 and L ∞ errors and the accuracies in the solution variables Dx, Dy, and Bz when we used a sequence of uniform meshes with resolutions ranging from 360 × 180, 720 × 360, and 1 440 × 720 zones. The results were compared with the suitably down-sampled results from the 2 880 × 1 440 zone mesh to extract the accuracy. We also ran a sequence of simulations on a two level composite mesh, where the global mesh as well as the refined patch had resolutions of 180 × 90, 360 × 180, and 720 × 360 zones. In other words, these runs had an effective resolution of 360 × 180, 720 × 360, and 1 440 × 720 zones, same as the uniform mesh runs. We see that our composite mesh with an effective resolution of 1 440 × 720 zones (i.e., with 720 × 360 zones on the global mesh and 720 × 360 on the refined patch) has errors that are roughly comparable to our 1 440 × 720 zone uniform mesh. Furthermore, second order of accuracy has been clearly reached on both sequence of runs. In the table, we also show the undivided divergence of the electric displacement, divided by the RMS value of the electric displacement. Notice that it is well-contained and close to machine zero in all instances. Notice that the finer meshes take many more timesteps than the coarser meshes, so there is a tendency for slightly greater divergence build-up on finer meshes than on coarser meshes, but again, that build-up is well-contained. Figures 13a-c show the pointwise error (relative to the maximum value of that variable) in Dx, Dy, and Bz when a composite mesh solution (720 × 360 zone coarse and fine meshes) is compared to a down-sampled uniform mesh solution (2 880 × 1 440 zones) using a second-order scheme. Figure 13d also shows the undivided divergence of the electric displacement relative to the RMS value of the same. We see that the divergence has been preserved to machine accuracy on the composite mesh. The final time of the simulation is shown. Table 14 is analogous to Table 13, but it has been carried out for the third order accurate runs. The same sequence of meshes was used. We see that third order of accuracy has been reached. Figures 14a-c show the pointwise error (relative to the maximum Fig. 11 Variation of the relative permittivity in the two-dimensional domain spanning [-14, 14] × [-7, 7]. It shows a cylinder with enhanced permittivity which is impacted from the left by a Gaussian pulse of electromagnetic radiation. When the mesh refinement is used, we use a refined patch that is centered on the cylinder and extends from [-7, 7] × [-3.5, 3.5], again with square zones that have linear dimension that is half that of the coarse mesh. The solid lines inside Fig. 11 show the extent of the refined patch value of that variable) in Dx, Dy, and Bz when a composite mesh solution (720 × 360 zone coarse and fine meshes) is compared to a down-sampled uniform mesh solution (2 880 × 1 440 zones) using a third order scheme. Figure 14d also shows the undivided divergence of the electric displacement relative to the RMS value of the same. We see that the divergence has been preserved to machine accuracy on the composite mesh. The final time of the simulation is shown. Table 15 is analogous to Table 13, but it has been carried out for the fourth order accurate runs. At fourth order, the accuracy is attained at 720 × 360 zones. Because the 2 440 × 1 440 zone simulation is extremely long-running, at fourth order, we use the 1 440 × 720 zone simulation as a reference run. We see that fourth order of accuracy has been reached. Interestingly, we see that the fourth order scheme is so accurate that it reaches its design accuracy on coarser meshes than the second-and third-order schemes. This highlights the value of a higher order accurate scheme. Figures 15a-c show the pointwise error (relative to the maximum value of that variable) in Dx, Dy, and Bz when a composite mesh solution (360 × 180 zone coarse and fine meshes) is compared to a down-sampled uniform mesh solution (1 440 × 720 zones) using a fourth order scheme. Figure 15d also shows the undivided divergence of the electric displacement relative to the RMS value of the same. We see that the divergence has been preserved to machine accuracy on the composite mesh. The final time of the simulation is shown. Observe that these errors are much smaller for the fourth order runs than the errors in Fig. 13, even though the mesh used in Fig. 15 has half the resolution as the mesh in Fig. 13. Notice too that there is no mesh imprinting in the error plots at the fine-coarse interfaces. In other words, there is no excessive pile-up of error at the fine-coarse interfaces, which is what one would have if there were a late-time instability. Furthermore, notice that these results have been obtained by an application of our standard algorithm over thousands of timesteps as encapsulated in Sect. 7 and Figs. 8, 9, and 10. There was no artificial "fix" applied to the algorithm like the fixes that previous practitioners have tried for curing the late-time instability in CED. It shows that we have obtained a conceptually tight solution of this long-standing problem.

Conclusions
In this paper we have addressed, and found, two very innovative resolutions of the problem of accurately prolonging a divergence-constrained vector field across successively refined meshes in an AMR hierarchy. While second-order accurate solutions of this problem were first invented in Balsara [1], the present work goes well beyond second order of accuracy. PDEs, where this problem arises include CED and MHD. In all such applications, the components of the vector field are collocated at the faces of a zone. This face-centered collocation is needed, because the update strategy for the numerical scheme requires such a collocation to mimetically reflect the constraints in the PDE. In subsequent papers, we will show how this prolongation strategy integrates with end-to-end adaptive solution of such PDE systems.
The problem of prolongation of constrained vector fields is indeed shown to be very challenging as the desired order of accuracy is increased. At up to fourth order of accuracy, an analytically exact solution has been found in Sect. 2 and the mathematical detail has been made explicit in Supplement B. However, with increasing order, this process becomes progressively more challenging, because we have to match not just a larger number of  Fig. 11. We used a second-order scheme with a 2 880 × 1 440 zone mesh. At that resolution, any further change in the solution becomes indistinguishable by eye. The solution on this mesh was used as a reference solution for all our second-order runs. Similar third and fourth order simulations were used as reference simulations for our third and fourth order tests. a-c Solution variables at time 0. d-f Solution variables at 3.5 × 10 -8 s, at which time the wave is interacting strongly with the refractive disk. g-i Show the solution variables at 5.7 × 10 -8 s, by which time the Gaussian pulse has scattered off the refracting cylinder. The black circle shows the extent of the disk with a radius of unity modes at the boundaries of the mesh but we also have to match a substantially larger number of constraints. The way out of this impasse was found in the almost divergence-preserving algorithm-WENO-ADP. The algorithm is made simple by the fact that one can match the modes in the faces in an order-by-order fashion. The result is a reconstructed vector field that spans the volume of a zone and, nevertheless, matches the moments that reside on the faces of the mesh. This reconstructed vector field is not exactly constraint-preserving. However, we have shown that it is almost constraint-preserving. In other words, the errors in the divergence constraint will still have an order property and the violation of the constraint goes down with increasing refinement in an orderly fashion. This suggests that an easy touch-up procedure can be invented that restores the divergence constraint exactly. This realization gives rise to the WENO-ADPT algorithm.
The WENO-ADPT algorithm is very versatile. First, it is very easy to extend to increasingly high orders and we have made the math explicit at up to fifth order. Second, our algorithm is very easy to implement. In fact, Sect. 3 shows that it is no more difficult to implement than a one-dimensional WENO scheme! Third, the touch-up procedure in Sect. 4 is provided in all its details. As a result, for prolongation of divergence-constrained vector fields on AMR meshes with a refinement ratio of two-which is the most popular caseone has just to transcribe the formulae provided here into code. Fourth, for refinement ratios greater than two, a very standard CLSQ procedure can provide the touch-up, and it too has been shown to work well.
In the course of doing this work, we have also been able to discover very efficient ways of implementing finite volume WENO schemes on structured meshes. This has been described in Supplement A, which is also foundational to many of the other sections in this work. In particular, we have been able to show that smoothness indicators for finite volume WENO schemes can also be written down as the sum of perfect squares. This advance mirrors a similar advance from Balsara et al. [6] as it pertained to finite difference WENO. As   Fig. 13 a-c Pointwise error (relative to the maximum value of that variable) in Dx, Dy, and Bz when a composite mesh solution (720 × 360 zone coarse and fine meshes) is compared to a down-sampled uniform mesh solution (2 880 × 1 440 zones) using a second-order scheme. d Also shows the undivided divergence of the electric displacement relative to the RMS value of the same. We see that the divergence has been preserved to machine accuracy on the composite mesh. The final time of the simulation is shown a result, even finite volume WENO methods have been made more robust and efficient with the innovations presented here. We also show that such finite volume WENO methods can be used to make very high order prolongation of scalar variables on mesh hierarchies. The full divergence-preserving algorithm has been presented in Sect. 7. We have also constructed a detailed test problem in Sect. 9. Via numerical tests we have shown that the methods presented naturally resolve the late-time instability that is well-known in the field of CED. Our results are order preserving and meet their design accuracies on composite meshes, as shown from Tables 13, 14, and 15. Furthermore, the divergence has been held down to machine accuracy on composite meshes, as seen from Figs. 13, 14, and 15. It is interesting to recall that careful attention to the second-order constrained prolongation problem in Balsara [1] led to novel methods for MHD (Balsara [2]) and CED (Balsara et al. [11]). In later papers we will also explore the algorithmic implications of the methods developed in this paper. Fig. 14 a-c Pointwise error (relative to the maximum value of that variable) in Dx, Dy, and Bz when a composite mesh solution (720 × 360 zone coarse and fine meshes) is compared to a down-sampled uniform mesh solution (2 880 × 1 440 zones) using a third order scheme. d Also shows the undivided divergence of the electric displacement relative to the RMS value of the same. We see that the divergence has been preserved to machine accuracy on the composite mesh. The final time of the simulation is shown