A Universal Predictor‐Corrector Approach for Minimizing Artifacts Due To Mesh Refinement

With nested grids or related approaches, it is known that numerical artifacts can be generated at the interface of mesh refinement. Most of the existing methods of minimizing these artifacts are either problem‐dependent or numerical methods‐dependent. In this paper, we propose a universal predictor‐corrector approach to minimize these artifacts. By its construction, the approach can be applied to a wide class of models and numerical methods without modifying the existing methods but instead incorporating an additional step. The idea is to use an additional grid setup with a refinement interface at a different location, and then to correct the predicted state near the refinement interface by using information from the other grid setup. We give some analysis for our method in the setting of a one‐dimensional advection equation, showing that the key to the success of the method depends on an optimized way of choosing the weight functions, which determine the strength of the corrector at a certain location. Furthermore, the method is also tested in more general settings by numerical experiments, including shallow water equations, multi‐dimensional problems, and a variety of underlying numerical methods including finite difference/finite volume and spectral element. Numerical tests suggest the effectiveness of the method on reducing numerical artifacts due to mesh refinement.


10.1029/2023MS003688
2 of 23 Despite the various choices of meshes and numerical methods, a commonly faced challenge for both the two types is the numerical artifacts or spurious reflections encountered at the place where there is a sharp change of mesh resolution.This issue becomes especially prominent for the non-conforming case, in which abrupt resolution changes are common by non-conforming faces or hanging-nodes.
The aim of this paper is to propose a method that can (a) minimize the artifacts due to mesh refinement and also (b) have a wide range of applicability to most of the numerical methods and underlying PDEs.To put our work in perspective, we review some of the existing work on numerical artifacts by mesh refinement.Early work dates back to at least the 1970s (e.g., Browning et al., 1973;Ciment, 1971Ciment, , 1972;;Phillips & Shukla, 1973) and to the 1980s (e.g., Berger, 1985Berger, , 1987)), including a series of work of Vichnevetsky (1981Vichnevetsky ( , 1985Vichnevetsky ( , 1987)), where wave propagation and spurious reflections were investigated using time-Fourier analysis.This issue was also covered in the review article of Trefethen (1982) on group velocity and the work of Celep and Bažant (1983) on elastic waves.Later work includes the study by Ullrich and Jablonowski (2011), where the decay rate of spurious reflections was investigated under different numerical settings.In brief, various techniques exist for minimizing artifacts by mesh refinement.However, these techniques are often tailored to specific scenarios and can be hard to generalize to other cases.For instance, upwind schemes are free from artifacts for the specific case of the advection equation, thanks to the use of only upstream information, but they still suffer from artifacts for general hyperbolic systems (Ullrich & Jablonowski, 2011).Artificial dissipation may remove spurious oscillations but can also hamper the conservation property or lead to a more stringent CFL stability condition.In addition, for low viscosity flow, the solutions are often sensitive to artificial dissipation, restricting the possibility of adjusting it only for the sake of minimizing artifacts.
Motivated by the previous discussion, in this paper we propose a simple predictor-corrector method to minimize artifacts by mesh refinement.The method introduces an additional grid setting with a different refinement interface, and then uses the partner simulation on this additional grid as a corrector for the original dynamics.It can also be related to data assimilation approaches (Chen & Stechmann, 2019;Torchinsky & Stechmann, 2023) in the sense that the partner simulation can be viewed as an "observation."This general framework allows the method to be easily applied, by construction, to most of the existing numerical methods, as an extension to their existing codes.Since the method does not require any modification of the existing methods, but simply serves as a corrector/post-processing at each time-step, this property makes it especially convenient to be used for the codes for which several options of numerical methods or physics parameterizations can be turned on/off or used in different combinations.For instance, this is the case for many weather/climate prediction models such as the Community Atmosphere Model (Dennis et al., 2012;Thomas & Loft, 2005), the Nonhydrostatic Unified Model of the Atmosphere (Marras et al., 2016), and the Eulerian/Lagrangian model (Prusa et al., 2008).In addition, due to the chaotic and turbulent nature of the weather/climate system, ensemble methods are frequently used in these models with multiple simulations running simultaneously, which could provide a setting in which our method could be applied to complex systems.While the use of a (full) partner simulation can be a valuable option for its simplicity, the doubling of the computational cost may be undesirable for some applications, and an alternative setup may more appropriate.As an alternative setup, one can use a partner mesh only over a small portion of the domain that is localized near the refinement interface, similar in spirit to the localization of the level-set method near an interface (Peng et al., 1999) or to the transition or buffer zone used in some grid nesting approaches (Davies, 1976;Davies & Turner, 1977;Debreu & Blayo, 2008;Giorgi et al., 1993;Lo et al., 2008;Marbaix et al., 2003;Phillips & Shukla, 1973;Urrego-Blanco et al., 2016), so the additional cost could be smaller.
A comparison of some different approaches and different mesh refinement configurations is shown in the schematic diagrams of Figure 1.One approach is to use a fine grid spacing in one region and a coarser grid spacing in another region (Figure 1) (e.g., Trefethen, 1982;Ullrich & Jablonowski, 2011).Another approach is to allow an overlapping region or feedback zone where the fine and coarse meshes overlap with each other (Figure 1b), and within which the fine-and coarse-mesh solutions can be nudged or relaxed toward each other (e.g., Davies & Turner, 1977;Marbaix et al., 2003;Urrego-Blanco et al., 2016).The overlapping region allows for a more gradual transition in the solution between the fine-mesh and coarse-mesh regions.In the predictor-corrector approach that is proposed herein (Figure 1c), the main goal is to correct the solution near a location of mesh refinement.To do so, a partner mesh is used and is configured to have a uniformly spaced mesh where needed.The partner mesh can be the same as the original mesh, except with a deliberate displacement of the location of mesh refinement.In this way, every spatial location in the domain will be covered by at least one mesh that has a uniform grid spacing in its local vicinity.

Numerical Artifacts and the Corrector
In this section, we briefly review the cause of numerical artifacts/spurious reflections and then propose the predictor-corrector method in the setting of a one-dimensional advection problem.The simplicity of this model problem allows us to more clearly demonstrate our method and explore its properties.The key idea of the method can be easily generalized to other problems and high-dimensional cases, which we shall explore in Section 4 and Section 5, respectively.

The Cause of Artifacts
Let us begin by considering a linear advection on one-dimensional grids.The equation reads Figure 1.Comparison of schematic diagrams for different approaches to mesh refinement and mesh configurations.(a) Simple setup of a mesh with fine grid spacing in one region and coarser grid spacing in another region (e.g., Trefethen, 1982;Ullrich & Jablonowski, 2011).(b) Transition layer with a fine mesh and a coarse mesh that overlap in a transition region, in which the solutions on the two meshes can be nudged toward each other (as indicated by vertical arrows) (e.g., Davies & Turner, 1977;Marbaix et al., 2003;Urrego-Blanco et al., 2016).(c) The predictor-corrector approach that is proposed herein, which like panel (b) involves an overlapping region, but where the solution at a regularly spaced mesh location is used to correct the solution near a location of mesh refinement.Also see schematic diagrams below in Figures 3 and 14.

10.1029/2023MS003688
4 of 23 where a > 0 is a constant representing an underlying right-going velocity.
To obtain a numerical solution of 1, the first step is to discretize the space: . . .< −2 < −1 < 0 < 1 < 2 < . . ., and we denote I j = [x j−1 , x j ] and Δx j = |I j |.The situation we are interested in is as follows: Namely, a sharp change in the grid spacing occurs at x 0 , where we may observe numerical artifacts.Based on the grid, a semi-discretization in space can be obtained.Here, we consider the FV methods, for which we have where u j represents an average of u on I j and F j represents an approximation of the flux at x j .The properties of the method are determined by how we choose the numerical flux F j .By choosing the upwind type flux F j = aF(…, u j−1 , u j ) we can suppress the artifacts since only upstream information is used.However, this approach is limited to this model problem.For more general problems (e.g., compressible gas dynamics, shallow water equations), both directions of information have to be used, and for this and other reasons such as ease of implementation, it is sometimes more convenient to make the choice of a central-type flux F j = aF(…, u j−1 , u j , u j+1 , …) (e.g., Nessyahu & Tadmor, 1990).
Unfortunately, this type of schemes allows spurious solution traveling in the wrong direction (left-going) that is not allowed in the PDE itself.This can be easily seen by investigating the dispersion relation of 2. Assuming the flux and considering locations away from x 0 , and by substituting the u j in 2 with   (−  ) , we have where h is the mesh size, ω is the frequency, κ is the wavenumber, and  d d represents the group velocity which governs the energy propagation (Trefethen, 1982) ) | .Note that the left and right-going waves are decoupled on the uniform grid such as the one used in Figure 2, but they would be coupled at the place of the refinement/coarsening if refinement/coarsening were used.Consequently, an artifact traveling leftward would be generated at any place where a sharp change of resolution occurs (Vichnevetsky, 1981).

Two-Grids Predictor-Corrector
Now we introduce the key idea of the predictor-corrector method.In the previous subsection we have shown that the artifacts are generated at the place of a sharp change of grid resolution.Therefore, if we have two grids with different refinement locations, then each of them can be used as a corrector for the other partner simulation.To be more specific, suppose we have two grids: satisfying x j − x j−1 = h r for j > 0 and x j − x j−1 = h l for j ≤ 0, and    −  −1 =  ℎ for j > 0 and    −  −1 =  ℎ for j ≤ 0. By letting  |0 − x0| =   0 , or in other words, by deliberately mismatching the resolution-changing positions of the two grids, we can create two information exchange channels with radius about R, such that in these channels, the correct dynamics on uniform grids can be used as correctors for the spurious dynamics.
Consider the case that h l < h r and  h < h .Suppose a wave packet starts from the very left end and travels rightward on both grids with the same initial setting.In this case, the wave shall first reach x 0 , at which it is crossing from the fine grid to the coarse grid on the mesh {x j }, while on the mesh  { x} , the wave packet is fully supported by a uniform fine grid.This corresponds to the corrector_1 in Figure 3.We can categorize this corrector as type DU AND STECHMANN 10.1029/2023MS003688 5 of 23 f2c-bf, which represents traveling from fine(f) grid to 2 coarse(c) grid using the background(b) uniform fine(f) grid information.Similarly, when the wave reaches   0 , the corrector_2 will be used and it has the type of f2c-bc.Since we have considered an advection problem with a specific grid setting (h l < h r and  h < h ), these two types of correctors (type f2c-bf and f2c-bc) are sufficient.However, for more general problems where waves can travel in both directions or for other grid settings, we have to also consider the other two types of correctors, namely, c2f-bc and c2f-bf.Their meanings are self-explanatory.Let us remark that for general wave problems, since waves can travel in both directions, the corrector types c2f-bf and f2c-bf should be regarded as a single type, and similarly for c2f-bc and f2c-bc.However, in the setting of advection, we can decouple them into four types, which allows us to examine them in more detail for better understanding.Now we introduce how to choose the correctors specifically.After choosing a time-stepping method (explicit or implicit) for 2, we obtain a full-discretization scheme: Behavior of wave packets of different wave-lengths: 8, 4, and 2.67 hr.In these simulations of 2, the initial condition is a wave packet at the center x = 1, defined by exp(−100(x − 1) 2 ) cos(κ(x − 1)) with κ = 50π, 100π, 150π, and the solution is plotted after evolving to time T = 1.The simulations use a central-flux finite volume method with fourth-order Runge-Kutta (RK4) time-stepping, and a uniform grid spacing.The parameters used are a = 1, h = 0.005, and dt = 0.5 hr., where we denote  () the polynomial space of degree k on the interval I. Furthermore, let P V and  P Ṽ denote the L 2 projections onto these two solution spaces, respectively.As an example, if ,1] .Given two weight functions   Ã ∈  2 (ℝ) taking values in [0,1], we introduce the corrector step as follows: These weight functions determine how much "trust" will be put into the partner simulations.The support of these functions should be limited by the distance between x 0 and   0 , since otherwise incorrect dynamics would be used as correctors.For instance, see Figure 3, if the radius of the support for the corrector_2 is larger than   0 − 0 , then the artifacts generated at x 0 on the top grid will be added to the solution on the bottom grid.

Weight Functions
In this section, we will explore the relations between the choices of the weight functions and their effects on minimizing the artifacts.The section is divided into three parts.In the first subsection, we explain the mechanism of artifact-minimizing by the corrector and give a formula (Equation 6) to estimate the correcting ratio, which is defined to be the ratio between the magnitudes of the artifacts with and without using the corrector.The analysis is carried out for the corrector type f2c-bf/c2f-bc since they allow us to exclude the artifacts caused by other sources.In the second subsection, we show that new types of artifacts can be generated by a sharp change of the weight functions.A formula is given to estimate the size of these artifacts (Equation 8).In the last subsection, we combine the results of the previous two parts and conclude on the guidelines of choosing the weight functions.

The Mechanism of Artifact-Minimizing
Let us start by investigating the corrector type f2c-bf.Namely, the corrector corrects the spurious artifacts for the case of a wave traveling from fine to coarse grid by using the background uniform fine grid information.
For simplicity, for now we assume the weight function w ≡ w m is a constant function in a left-neighborhood (with length L cr ) of x 0 , at which the grid suddenly changes its resolution.To help the explanation, we use Figure 4 for visualization.When a right-going wave is crossing x 0 , for each time-step, artifacts with wavelength lying in (2h, 2h + ϵ) are generated.By 3, these artifacts, which are represented by the oscillation in Figure 4, have a negative group velocity.Therefore, we can decompose the numerical solution into the right-going wave and the left-going artifacts.On the other hand, on the uniform fine mesh, the wave is crossing x 0 without any artifacts being generated.The above concludes the predictor step.For the corrector step, both the right-going wave and the artifact are multiplied by the weight function (1 − w m ).On the background grid, the right-going wave is multiplied by w m and added to the solution.This step corrects the missing w m portion of the numerical solution caused by multiplying (1 − w m ).Most importantly, the artifact has been multiplied by (1 − w m ).This suggests that the artifact should decay exponentially as it travel leftward until it escapes the support of the weight function.Based on this, we can have an estimate for the ratio between the sizes of the artifacts with and without using the corrector: where v s represents the group velocity of the artifact.We will call this constant the correcting ratio.Similar arguments can be made for the corrector type c2f-bc.
We next present an experiment to verify 6.The grid nesting is set as follows.For the two grids, we set   = 4 = Ã = Ã = 0.005 .Namely, for the first grid, the size of the coarse grid on the right is four times of the fine grid on the left; the second grid is chosen as a uniform fine grid and will be used for the corrector.For the solution, we choose a = 1, the exact solution to be u(x, t) = sin(8πx − t), and use an oscillating Dirichlet boundary at x = 0 and an absorbing boundary at x = 2.We shall take snapshots at t = 2 when the solution can be observed to have obviously entered the time-harmonic region.
DU AND STECHMANN 10.1029/2023MS003688 7 of 23 We carried out four tests, where we vary the three determining factors of r cr that appear in 6, namely, the magnitude and the support length of the weight function w m and L cr , and the length of the time-step Δt; these values are collected in Table 1.For all simulations, we set  w = 0 so the simulation on the uniform fine grid is independent of its partner simulation.The results of these simulations are collected in Figure 5.Note that we have chosen very small values for w m such that the artifacts would not be too small to be observed.

Artifacts Caused by Weight Functions
In the previous subsection, we have shown that the artifacts generated by mesh refinement can be efficiently reduced by the corrector in an exponential order; see 6.However, as we will show in this subsection, the weight function itself can generate other artifacts.See Figure 6, for which we have used the same experimental setting for Figure 5, except now we also include the test for the corrector type f2c-bc to compare.In the left column, which is for the corrector type f2c-bf, we observe that the artifacts decay quickly thanks to the relative large value of the weight function w.On the other hand, for the corrector type f2c-bc, we observe artifacts being generated at the left endpoint of the support of the weight function.
The artifacts shown in Figure 6, right column, are caused by the sharp change of the weight function.We next derive a formula to estimate the amplitude of these artifacts.First note that we can rewrite the corrector Equation 5 as follows:      agrees well with u, no correction is applied in the neighbor of x 0 − L cr so the artifacts can be negligible (this corresponds to the corrector types f2c-bf and c2f-bc; see the left column of Figure 6).However, when    has a relatively large difference from u, this added quantity  P (( ũ − )) has to be taken into account, which corresponds to the corrector types f2c-bc or c2f-bf.In what follows we focus on these cases.Note that the artifacts with negative group velocity are generated in a neighbor of x 0 − L cr by the jump of the weight function from 0 to w m .Suppose all the added jump  P (( ũ − )) close to x 0 − L cr contributes to these artifacts, and we assume an energy balance between the added energy and the energy traveling leftward, then we can have an approximate estimate for the amplitude of the artifacts generated at x 0 − L cr : where v sw is the group velocity of the artifacts generated at x 0 − L cr , and     is a constant that depends on the mismatch between the solution u and    .We next carry out an experiment to verify 8.We use the same experimental setting as the one used for Figure 5, except now we use a coarse grid corrector; namely, we choose  4ℎ = ℎ = h = h .We vary the value of w m and the size of the time-steps and then observe the changes of the amplitude of the artifacts.See Figure 8; the experiment results agree with the estimate by 8. Namely, the artifacts caused by the weight functions are proportional to w m but anti-proportional to Δt.

Guidelines for Choosing the Weight Functions
For general corrector types, we need to consider the artifacts caused by both the mesh refinement and the weight functions.Combing Equations 6 and 8, we have an estimate for the final outcome of the artifacts: where C 1 and C 2 are constants independent of w m and L cr .Group velocity of the artifacts are assumed to be 1 for simplicity.By Equation 9, we know that increasing the value of L cr can decrease the size of the artifacts.Note that L cr should principally be limited by the distance between the two refinement interfaces (see the discussion at the end of Section 2.2).This fact would motivate us to choose a larger distance between the two refinements to have a better effect on minimizing the artifacts.Though, let us remark that a larger support can increase the communication burden between the simulations and hence the support length L cr needs to be chosen accordingly based on the demand.On the other hand, since the first and the second term in 9 are negatively and positively correlated with w m , respectively, the value of w m has to be chosen in a way that minimize the contribution from both the two terms.
Considering that the second term of 9 is caused by a sharp change of the weight function, namely the jump of w from 0 to w m , an alternative way of choosing the weight that can possibly decrease the error from this term is to consider a smooth transition of the weight: We keep the maximal value w m and the support length of the weight functions L cr unchanged so that we can fairly compare the results.By comparing Figure 9 and Figure 5, we observe similar decaying pattern of the artifacts by using the tent-shape weight function given by 10.Though, the decaying is less obvious.This is expected since the value of 1 − w becomes closer to 1 as the artifacts travel leftward and approaches x 0 − L cr , leading to a weaker effect of the corrector (locally adopting 6).On the other hand, by comparing Figure 10 with Figure 8, one can see that the errors from the artifacts have been reduced by roughly a factor of 20.Hence, we observe that the artifacts caused by weight functions have been significantly decreased, thanks to the smooth transition of the weight at x 0 − L cr .
Not surprisingly, the tent-shape weight function generates much smaller artifacts thanks to its relatively smooth transition from 0 to the maximal value w m .Despite that it is less effective on suppressing the artifacts caused by mesh refinement, we find that we can always choose a larger value of w m so that the overall effect of minimizing the artifacts is better than the case of using a step-function as the weight.While using the step-function as weight function may not provide the best results in practice, it is useful for providing theoretical guidance and understanding.
To conclude, we find that using a tent-shape function as the weight function achieves better results than using a step-function on the overall effect of minimizing the artifacts, at least for the corrector types f2c-bc and c2f-bf, in which cases the artifacts introduced by the weight functions themselves become the prominent issue.In addition, increasing the support length of the weight L cr can always decrease the artifacts caused by mesh refinement.But it needs to be limited by the distance between the two refinement locations of the two grids, and also by the demand to save communication costs between the two simulations.Finally, increasing the maximal value of the weight function w m can decrease the artifacts due to mesh refinement, but can also increase the additional artifacts by the weight function itself.Therefore, w m needs to be optimized to achieve the best effect.For tent-shape weight functions, our experiments suggest that choosing the slope of the function     cr to be around 0.5 can achieve a good balance and provide near optimal results of minimizing the overall artifacts, provided that L cr is not too small.
At the end let us remark on the other possible choices of the weight functions, such as trigonometric-type functions and higher-order polynomials.Our preliminary experiments suggest a similar artifacts-minimizing effect compared to the linear-type weight represented by 10.Consequently, we do not delve into the more details about them in this paper, due to their seemingly inconsequential nature and the consideration of the length limitation.

Tests for Shallow Water Equations
In this section, we move beyond the advection equation and present numerical experiments to test the predictor-corrector method for shallow water equations.These equations allow waves traveling in both directions and contain nonlinear terms, and consequently represent a more realistic setting for applications.In addition, some of the traditional techniques of removing artifacts by mesh refinement fail to work for shallow water 10.1029/2023MS003688 12 of 23 equations (see Ullrich and Jablonowski (2011)), which further motivates our investigation of applying the predictor-corrector method to these equations.We begin by considering on the domain [a, b] with periodic boundary conditions.Here h(x, t) is the pool depth/height, m(x, t) is the momentum, and g is a parameter that represents the gravitational acceleration.We consider the discretization of and use RK4 for time-stepping with the length of the time-step Δt.Again, we have used u j to represent the average of the solution on the interval [x j−1 , x j ], and used F j to represent the flux at the location x j .The above concludes the predictor step.For the corrector step, the weight functions w and  w are chosen as the tent functions (dashed lines) shown in Figure 11.We assume all the tent functions have the same radius of support L cr and the maximal value w m .These notation first appeared in 10.We next carry out the numerical tests.The initial solution is chosen to be a Gaussian-type perturbation at the center of the domain [a, b] based on the reference height  ℎ ≡ 2 : We specify all the aforementioned parameters in this section using the values collected in Table 2.
In the first test, we set the fine grid in the interior and the coarse grid in the exterior regions.In this setting, the initial Gaussian distribution at the center will first decompose into a left going and a right going distribution, and then leave the interior fine resolution grid and enter the exterior regions of coarse resolution.We compare the results without and with using the corrector and collect them in Figure 12.We observe that the numerical artifacts have been mostly removed.We remark that in this simulation, the corrector type f2c-bf was used at the location m r = 3 on the top grid {x j } and    = −1 on the bottom grid; the corrector type f2c-bc was used at the location    = 4 and m l = −2.
In the second test, we swap the resolution in the interior and the exterior regions of the first test.In this case, the solution will travel from the coarse grids to the fine grids, and the correctors with the types c2f-bc and c2f-bf will be triggered.The results are collected in Figure 13.We again observe that the corrector has effectively removed most of the artifacts.Therefore, the method appears to work well in a case where nonlinear waves can simultaneously propagate in multiple directions, both leftward and rightward, with all of the four types of the correctors being used.

Mesh Setting and the Approximation Spaces
In this section, we introduce the predictor-corrector method for multi-dimensional problems.As a demonstration, we consider the implementation of the method for SE methods on 2D meshes.Note that 2D meshes allow hanging nodes which do not occur in 1D; see Figure 14.We consider two meshes of a rectangular domain Ω = [x 0 , x 1 ] × [z 0 , z 1 ] which is contained in the x-z plane and we denote them as   and   .Both meshes are composed of the fine meshes in the bottom regions and the coarse meshes in the upper regions, where the divisions between fine/coarse regions are given by the horizontal line z = z 0 for the mesh   and the line  = Ã0 for the mesh   .Again, we have chosen the setting such that  0 ≠ Ã0 so that we are able to use the dynamics on a region of uniform meshes to correct the spurious dynamics at the place where mesh refinement/coarsening or hanging nodes occur.

With the meshes
 and   , the approximation spaces by SE are defined as follows: where  { * ,, } =1,. . .,(+1) 2 are the coefficients of the numerical solution restricted on the element K, and  { } =1. . .(+1) 2 are the SE basis functions on K. Suppose   is the push-forward map from the reference element [−1,1] 2 to the physical element K, then we have , where  {} =1,. . .,(+1) 2 are the pre-determined SE basis functions on the reference

Implementation of the Predictor-Corrector
Now we explain how the predictor-corrector can be implemented for SE methods.Like the 1D case (see 4), the predictor step is determined by how we choose the spatial and temporal discretization.For instance, for spatial discretization, conforming methods like the traditional FE methods can be used, or we can use DG methods together with a choice of the numerical flux; for temporal discretization, Runge-Kutta (RK) or BDF methods can be used, either in their explicit or implicit forms.In all cases, the predictor step can be abstracted as the following updating algorithm: where    * ,ℎ with * = a or b represents the numerical solution at the time-step t = nΔt.To be more specific,

𝑗𝑗=1
* ,, , (, ) approximates the exact solution u(nΔt, x, z).Then the corrector step is defined as follows: The implementation of the above corrector involves the calculation of certain weighted L 2 projections at the two narrow bands around z = z 0 and  = Ã0 , which are the support regions of the weight functions w a and w b ; see the gray regions in Figure 14.Equation 13a corresponds to the leftward arrow and the gray region in the left mesh in Figure 14.It involves the calculation of the element-wise weighted L 2 projections with the weight w a and 1 − w a , and also the calculation of a weighted-gathering projection since it takes the degrees of freedom (DOFs) of four elements and combine them to update the DOFs on one element.This weighted-gathering projection corresponds to the upper part of the gray region on the left.On the other hand, Equation 13b is represented by the rightward arrow and the gray region on the right.It involves the calculation of a new type of weighted-restricting projection operator, which corresponds to the bottom part of the gray region on the right.In conclusion, there are three types of projections that need to be calculated at the narrow bands near z = z 0 and  = Ã0 : (a) element-wise weighted projections, (b) weighted gathering projections, and (c) weighted restricting projections.
To explain how these three types of projections can be calculated for SE methods, we first introduce some notation.Let  {   ,  } (+1) 2

𝑗𝑗=1
be the quadrature points and weights on the reference element [−1,1] 2 .We use λ to represent the weight functions for correctors and it could be w a , w b , 1 − w a , or 1 − w b .The calculation of the element-wise weighted projections is the most straightforward.Suppose the element is K, whose four vertices, counting counter-clockwise from bottom-left to upper-left, are denoted as  )) . This flux can be regarded as the Rusanov flux with zero artificial viscosity, which allows us to exclude the source of artifacts dissipation by viscosity.This method can also be related to conforming methods using mortar elements techniques except that no averaging on vertices needs to be done.Here   nbr * ,ℎ is the projection of the solution of the neighbor elements of an element K to the boundary of K.This projection becomes a simple restriction when there are no hanging-nodes.If hanging-nodes are involved, we again need to calculate a gathering and a restriction projection operator, which can be regarded as the one-dimensional version of the gathering/restriction operators we have introduced in the previous subsections with a constant weight λ ≡ 1.We refer to Kopera and Giraldo (2014) for more on numerical techniques of dealing with hanging nodes in the SE setting.

On both meshes 𝐴𝐴
 and   , we initialize the numerical solution   0  and   0  by using the same exact solution of a Gaussian-type distribution 1+ 2 .The time-stepping method together with the spatial discretization given by 14 determine the predictor-step 12.
For the corrector-step 13, the weight function w a and w b are chosen as the tent-shape functions in z direction, with the center lines chosen as z = z 0 and  = Ã0 , and the radius of the support taken as L cr = 2; see Figure 17.
We next carry out two experiments.In the first experiment we choose μ = 0, which gives the advection in the perfect upward direction.See Figure 18; without the corrector, many high-frequency artifacts appear; see, for instance, the region with 2 < x < 4 and 0 < z < 15.When the corrector is used, on the other hand, the numerical artifacts have been removed to a degree that they can hardly be seen.To help see these artifacts, we present another visualization along the z-axis; see Figure 20.We observe that the artifacts have been significantly reduced to a tiny amount.Note that the support of the weight function only contain one big and two small elements (along the z-axis direction); see Figure 14.This suggests that our method has automatically taken advantages of the high order elements, in comparison to the low order FV cases in the previous sections, where often around 10 times more elements need to be included in the support to achieve a similar artifacts-minimizing effect.
In the second test, we increase the value of μ to be μ = 0.1, which gives an advection velocity such that the Gaussian-type distribution will hit the mesh refinement interface with a tilted angle; see Figure 19.In this case, without the corrector, one can again see many high-frequency artifacts; see, for instance, the region with 3 < x < 5 and 0 < z < 15.When the corrector is used, on the other hand, we again observe that the artifacts have been effectively removed.

Conclusion
In this paper, we have proposed a universal predictor-corrector method to minimize numerical artifacts caused by mesh refinement.Analysis is given to understand the mechanism of the method for a model advection problem.Numerical experiments are carried out to illustrate its general applicability, for different PDEs (e.g., the linear advection equation and nonlinear shallow water equations), for 1D and multi-dimensional problems, and for a variety of different numerical methods (e.g., finite difference/FV methods and spectral elements).The numerical results suggest the effectiveness of the method for removing numerical artifacts.
The method here is based on the communication between the simulations on two meshes with different refinement locations.Considering that this communication only needs to be done at a small neighborhood of the coarse/fine mesh interface, its computational cost is in general small compared to the numerical schemes that update the simulations.In addition, the communication can be done after several time-steps instead of every time-step, which can further save computational costs.However, reducing the communication frequency can change the effect of the corrector, as is suggested by our preliminary numerical tests not included in the paper.The tests show that a smaller communication frequency can reduce the artifacts caused by the weight function but may reduce the correcting ratio.Therefore, the magnitude of the weight w m needs to be accordingly adjusted larger to have a better effect on reducing the overall artifacts.More detailed investigation of this topic and some others, such as the combination with ensemble methods, the exploration of other weight functions, exceeds the scope of this paper but may constitute future work.
Compared to existing methods such as applying a smooth transition layer for the mesh-resolution change (Moeng et al., 2007;Tang et al., 2013), and/or adding hyper-viscosity (Guba et al., 2014;Tang et al., 2013) to damp out unphysical artifacts, the proposed method here is more flexible with the mesh-resolution transition (e.g., allowing either an abrupt or smooth change in resolution), and it introduces less artificial viscosity and so it may lead to a better energy conservation property.
The tests here were conducted with two grids-an original grid and a partner grid-as a setup that is simple and requires minimal changes to a computer code.However, the use of two grids would double the computation cost in comparison to a single grid, if computed in serial.Also, if two grids are used for turbulent and chaotic fluid motion, synchronization would be needed occasionally to prevent the two simulations from diverging from each other.To overcome the extra computational cost that could arise from two grids, a different approach is to use a single grid setting, accompanied with a partner grid that is local and only exists in the vicinity of the mesh-resolution transition boundary.Such an alternative approach is similar in spirit to the localization of the level-set method by introducing a partner grid that is localized near an interface (Peng et al., 1999) or to the transition or buffer zone used in some grid nesting approaches (Davies, 1976;Davies & Turner, 1977;Debreu & Blayo, 2008;Giorgi et al., 1993;Lo et al., 2008;Marbaix et al., 2003;Phillips & Shukla, 1973;Urrego-Blanco et al., 2016) so the additional cost could be smaller.
Future work could also include the investigation of the methods here in the presence of complex physics parameterizations, such as earth system modeling, with geophysical fluids that are generally turbulent.Physics parameterizations, such as parameterizations of moist convection and cloud microphysics, are in use in a variety of forms in both idealized and comprehensive numerical simulations (Bendall et al., 2020;Dias & Pauluis, 2009;Khouider & Majda, 2005;McIntyre et al., 2020;Tissaoui et al., 2022;Wetzel et al., 2020;Zarzycki et al., 2015).Such parameterizations were beyond the scope of the present paper, where here the focus was given to the formulation and justification of the methodology for minimizing refinement artifacts, and focus was given here to testing other aspects such as different equations in different spatial dimensions and with different classes of numerical methods.Given the complex and varied form of many physics parameterizations, it can be difficult to design specialized strategies for minimizing refinement artifacts for individual parameterization strategies, and it may be promising to investigate the universal strategy proposed here.

Figure 3 .
Figure 3. Mutual corrector through deliberate displacement of the resolution-changing positions.The two channels are not necessarily non-overlapping as is drawn.

Figure 4 .
Figure 4. Mechanism of artifact-minimizing for corrector type f2c-bf.See the text for a detailed description.The procedure is repeated for every time step.At each step, the amplitude of the solution is maintained, while the amplitude of the artifact (or oscillations) shrinks by A osc ← (1 − w m )A osc by applying the corrector.

Figure 5 .Figure 6 .
Figure 5.Comparison of the correcting ratios for the corrector type f2c-bf by using different weight functions and time-steps.The simulations use central-flux finite volume with RK4 time-stepping.

Figure 7 .
Figure 7. Balance between the added energy from the jump of the weight function and the energy traveling leftward.

Figure 8 .
Figure 8.Comparison of the artifacts generated by weight functions for the corrector type f2c-bc, using different values of the weight and different sizes of the time-steps.

Figure 9 .
Figure9.Illustration of the correcting ratios from using a weight function that is tent-shaped.Same setting as Figure5except using 10 for weight functions.
[a, b]  into two grids {x i } and  { x} as is demonstrated in Figure11.Each grid uses a fine resolution of h f at the center region and a coarse resolution of h c in the exterior regions, or vice versa.On both grids, we discretize 11 by FV methods with central flux  + ( − −1)∕Δ = 0 where  = 1 2 (() + (+1)),

Figure 10 .
Figure10.Illustration of improved results and smaller artifacts from using a weight function that is tent-shaped.Same setting as Figure8except using 10 for weight functions.

Figure 11 .
Figure 11.Mesh setting and the weight functions (dashed lines) for a periodic domain.
19422466, 2023, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023MS003688,Wiley Online Library on [29/08/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 2 .They are the tensor products of the Lagrange basis functions based on the Gauss-Lobatto quadrature points.For more information about SE methods, see, for instance,Marras et al. (2016).

Figure 12 .
Figure 12.Comparison of the solution (water depth/height h) for shallow water equations without (left) and with (right) the corrector.Fine and coarse resolution are used in the interior and the exterior regions, respectively.The top and the bottom row represent the simulation on the grid {x i } and  { x} , respectively.

Figure 13 .
Figure13.Same settings as those used for Figure12except for using coarse resolution in the interior and fine resolution in the exterior regions.

Figure 14 .
Figure 14.Mesh setting for 2D.The left and the right meshes represent two different discretizations of the same domain Ω but are plotted separately for the ease of visualization.

Figure 15 .
Figure 15.Weighted restriction and gathering projection operators.

Figure 16 .
Figure 16.Initial conditions of the function u(t, x, z) at time t = 0 for the numerical solution on the mesh   (left) and the mesh   (right) with polynomial degree k = 4.The left and right panels differ in the z location of the mesh refinement interface; it is at z = 18 in the left panel and z = 20 in the right panel.

Figure 17 .
Figure 17.The weight functions w a (left) and w b (right).

Figure 18 .
Figure 18.Comparison between the numerical solution at t = 12 without (top) and with (bottom) the corrector.The solution propagates upward, in a direction that is orthogonal to the mesh refinement interface.

Table 1
Parameter Values Used for Tests With Different Weight Functions and Time-Steps 19422466, 2023, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023MS003688,Wiley Online Library on [29/08/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Table 2
Parameter Values Used for Tests With the Shallow Water Equations