{'fraction_non_alphanumeric': 0.030392992563416713, 'fraction_numerical': 0.014740895329355398, 'mean_word_length': 4.397667777248929, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient and Accurate Adaptive Resolution for Weakly-Compressible SPH

In this paper we propose an accurate, and computationally efficient method for incorporating adaptive spatial resolution into weakly-compressible Smoothed Particle Hydrodynamics (SPH) schemes. Particles are adaptively split and merged in an accurate manner. Critically, the method ensures that the number of neighbors of each particle is optimal, leading to an efficient algorithm. A set of background particles is used to specify either geometry-based spatial resolution, where the resolution is a function of distance to a solid body, or solution-based adaptive resolution, where the resolution is a function of the computed solution. This allows us to simulate problems using particles having length variations of the order of 1:250 with much fewer particles than currently reported with other techniques. The method is designed to automatically adapt when any solid bodies move. The algorithms employed are fully parallel. We consider a suite of benchmark problems to demonstrate the accuracy of the approach. We then consider the classic problem of the flow past a circular cylinder at a range of Reynolds numbers and show that the proposed method produces accurate results with a significantly reduced number of particles. We provide an open source implementation and a fully reproducible manuscript.


Introduction
It would appear that a meshless particle method would be naturally suited for adaptive resolution. However, accurate adaptive resolution for Smoothed Particle Hydrodynamics (SPH) in the context of weakly-compressible and incompressible fluid flow is still a challenging area of current research [1].
In the context of incompressible and weakly-compressible fluid flow, there have been some valuable developments starting with the pioneering work of Feldman and Bonet [2] where the particles are adaptively split in an accurate manner. This work has been extended further to include particle merging by Vacondio et al. [3,4] and applied to the shallow water equations [5], soil simulation [6], fluid-structure interaction [7]. The method has been designed to be very accurate and a great deal of care is taken when splitting and merging particles. However, the accuracy comes at a significant cost since each coarse particle splits into 7 particles in two dimensions and around 14 in three dimensions. This leads to an enormous increase in the number of particles as the regions are refined. The particle de-refining method merges particles pairwise and it is argued [8] that the method is computationally expensive since the rate of splitting particles is significantly larger than the rate of merging. While the method is designed to be accurate, the resulting refined particles also employ a very large smoothing radius in comparison to what would be expected in a fixed particle size discretization of the problem with a similar number of particles. This poses significant additional performance limitations on the method. Moreover, the method relies on manual specification of the spatial regions where the adaptation is desired. This is inconvenient in general and especially when the bodies are moving.
Barcarolo et al. [9], Sun et al. [10], and Chiron et al. [8] refine each coarse particle (also called parent particle), in two dimensions, into 4 child particles but also retain the coarse particle. The parent particles are passively advected in the refined regions. This implies that each coarse particle effectively splits into five particles. This reduces the number of refined particles when compared with [3]. The significant advantage with this approach is that de-refining particles is simple to implement; the parent particles are re-activated and the child particles are removed. This approach has also been used for some impressive multi-resolution simulations using the δ + -SPH scheme [10,11]. Another significant advantage is that the smoothing length chosen is much smaller than the typical values chosen in the approaches of [3]. In order to handle the interactions between the child and parent particles, a particle property γ, is added to each particle. In the transition regions this value is between 0 and 1 whereas in regions with uniform particle smoothing length, the value is either 0 or 1. When the value is zero for a particle, the particle is effectively switched off and when it is one it is active. Intermediate values allow for the use of both the parent and child particles.
Chiron et al. [8] further refined this method by taking inspiration from traditional Adaptive Mesh Refinement techniques to create an Adaptive Particle Refinement (APR) procedure. In the intermediate regions where larger particles are refined into smaller particles, both the parent and child particles are retained and only particles of the same size interact and the properties are carefully interpolated between the parent and child particles. The difficulty with the approaches of Barcarolo et al. [9], Sun et al. [10], and Chiron et al. [8] is that coarse particles effectively split into five particles in each level of refinement. Furthermore, there are additional complications due to the special handling required for the parent and child particles either by the use of the γ parameter or by the use of prolongation and restriction operations in the APR method. It is also not entirely clear what would happen in high-strain fluid flows where the four child particles would drift significantly apart away from the parent particle.
Recently, another approach for dynamic particle splitting and merging has been proposed by Yang and Kong [12,13] and applied to multi-phase fluid simulations. This approach is similar to that employed by Vacondio et al. [3] but each coarse particle is only split into two child particles. The parent particle is removed. However, the placement of the child particles is done carefully along the perpendicular bisector of the line joining the parent particle to its nearest particle. This method will only work in two dimensions and no procedure for the three dimensional case is proposed. The advantage with this approach is that the particle refinement is much more gradual without a very large increase in the number of particles. Merging is done only between two particles and therefore there is no profusion of particles. The proposed method also elegantly handles gradual refinement of the resolution around an interface using a single parameter. This has been demonstrated for multi-phase problems [13]. It appears that no detailed study of the accuracy of the method has been performed. However, previous accuracy studies by [2] suggest that splitting particles into only two child particles would introduce significant error into the solution. Moreover, the method has only been demonstrated for two dimensional fluid flows. An alternative to the distance-based spatial adaptation is the recovery-based a posteriori error estimator [7], where the error in the SPH velocity gradient is measured and particles are adaptively refined in regions where the error exceeds a tolerance.
In the area of computer graphics, Desbrun and Cani [14] use splitting and merging operators in the SPH method and applied it to highly deformable substances. Adams et al. [15] use extended local feature size to adaptively refine the particles in the regions of geometric interest. Solenthaler and Gross [16] use two-scales, a lower resolution and a higher resolution, and couple the two with appropriate boundary conditions and feedback forces. However, these works are designed more for computer graphics applications and do not test the accuracy with any standard benchmark problems.
Spreng et al. [17] have proposed the use of the method of Vacondio et al. [3] for performing adaptive particle resolution for structural mechanics problems. They also propose a novel method to merge multiple particles by considering neighboring particles which are identified in two different ways. It is not clear if the proposed algorithm is parallel as the details of the implementation are not discussed. In a more recent work, Spreng et al. [18] have proposed the use of adaptive refinement to improve the discretization errors.
Recently, Sun et al. [19] have employed the adaptive particle refinement and de-refinement approach to strongly compressible, multi-phase fluid flows. However, the main focus of the adaptive particle refinement is to ensure a homogenous and isotropic particle distribution when the fluid is highly compressible.
The idea of splitting and merging particles is not new and has been successfully applied in the context of vortex methods [20]. This technique has also been used by various researchers employing SPH for computer graphics [14,15,16]. However, the challenge in implementing this with the SPH method for incompressible and weakly-compressible fluids is to have a method that is both accurate and computationally efficient with a minimum of numerical parameters. This is a significant challenge. It bears emphasis that none of the existing adaptive resolution schemes for fluid flow problems with widely varying scales [3,11,8] feature an automatic adaptation strategy, nor do they inherently support complex moving geometries or provide any ability to introduce solution-based adaptivity. The methods of [7,19,18] do support moving geometries and solution-adaptivity however, they do not seem to have been applied to problems with widely-varying scales.
In this paper we propose a new approach which is automatically adaptive, computationally efficient, accurate, supports moving bodies, and solution adaptivity. The basic strategy is to split and merge particles carefully as originally proposed by Feldman and Bonet [2] and Vacondio et al. [3]. However, we adaptively merge particles to reduce the large particle counts. This is done in a computationally efficient manner, in parallel, and our simulations suggest that this approach is also accurate. We are thus able to control the particle refinement adaptively so as to effectively only double the number of particles in each refinement region while retaining accuracy. In addition, we carefully set the smoothing radius of the refined particles to be optimal for the particular refinement region thereby further improving performance in comparison to the approach of Vacondio et al. [4]. We use ideas inspired from the work of Yang and Kong [13] to automatically set the refinement criterion. This allows us to specify the geometry, a few parameters determining the maximum and minimum length scales and the algorithm automatically refines the particles as required. We discuss in some detail the algorithm proposed and show how it can be used to (i) handle complex geometries, (ii) specify user-specified refinement regions, (iii) handle moving geometries, and (iv) be used for solution-based adaptivity. We do not extensively explore solutionbased adaptivity in this work but outline the basic ideas and demonstrate this with some simulations. The algorithms employed in this work are parallel and in principle may be executed on a General-Purpose Graphics Processing Unit (GPGPU).
We only consider two dimensional problems in this manuscript but in principle the ideas naturally extend to three-dimensional cases. Although we use a modified EDAC-SPH [21] scheme for the SPH discretization any similar method could be used. In the present work we do not consider any free-surface problems, however, our adaptive refinement algorithm can be easily extended to work with such problems. We consider several simple benchmark problems to demonstrate the accuracy of the approach. We then simulate the flow past a circular cylinder at a variety of Reynolds numbers in the range 40 -9500 and compare these with some very well established simulations to show that the method is capable of resolving the necessary details with a minimum of particles. This translates to a proportional reduction in the computational time. The new method allows us to perform such computations with far fewer particles than reported elsewhere with the SPH method. For the case of the flow past a circular cylinder the results we present require at least an order of magnitude fewer particles than those reported in [11] for a similar resolution. Finally, we note that none of the existing methods for adaptive SPH feature open source implementations. We provide a fully open source implementation based on the PySPH framework [22,23]. The source code can be obtained from https://gitlab.com/pypr/adaptive_sph. Our manuscript is fully reproducible and every figure is automatically generated through the use of an automation framework [24].

The SPH method
In this paper we deal specifically with weakly-compressible flows. We use the entropically damped artificial compressibility (EDAC) method [21] to simulate the weakly-compressible flows. The position update, pressure evolution, and momentum equations in the EDAC formulation are, where r, u, p, and t denotes the position, velocity, pressure, and time respectively. ρ is the density, ν is the kinematic viscosity of the fluid, c s is the artificial speed of sound, f is the external body force, and ν e is the EDAC viscosity parameter. In order to further enhance the uniformity of the particles we use the transport velocity formulation [25], with the corrections incorporated [26]. Then the above equations are re-formulated as, whereũ refers to the transport velocity, andd (.) dt = ∂(.) ∂t +ũ · grad(.) is the material time derivative of a particle advecting with the transport velocityũ. The computation of the transport velocity is shown in section 2.1.
Remark. In our numerical experiments with the Taylor-Green problem we found that the addition of the divergence correction terms in the pressure evolution equation is crucial for accuracy. However, we find that the use of the last two terms in the momentum equation (6) introduces noise where the particles are merged or split. Consequently, we do not use them in this work. We note that Sun et al. [27] observes that the effect of these terms in the momentum equation is minor.
We discretize the governing equations using variable-h SPH. The domain is discretized into points whose spatial location is denoted by r i , where the subscript i denotes the index of an arbitrary particle. The mass of the particle, which vary as a function of space, is denoted by m i , and its smoothing length by h i . In the variable-h SPH the density is approximated by the summation density equation using a gather formulation [28,5] written as, where, W (|r i − r j |, h i ) is the kernel function. We use the quintic spline kernel in all our simulations, the quintic spline kernel is given by, where σ 2 = 7/(478πh(r) 2 ), and q = |r|/h. The EDAC pressure evolution equation in variable-h SPH (see [29,5], for a derivation of the terms in the R.H.S) is given by, dp dt where ρ 0 is the reference density, p i is the pressure of particle i, ρ j is the density of the j th particle computed using summation density eq. (7), u ij = (u i − u j ), r ij = |r ij | = |r i − r j |, β i is the variable-h correction term [5], which in d dimensions is given by, P i and P j are given by, here we employ the pressure reduction technique proposed by Basa et al. [30], where, the average pressure is computed as, where N i is the number of neighbours for a particle with index i, and The EDAC viscosity of the pressure diffusion term in the EDAC equation with the SPH discretization is given by, where α e = 1.5 is used in all our simulations. Since this is a function of the smoothing length, which is varying in space, we use the approach of Cleary and Monaghan [31] to model the pressure diffusion term where, ν e,ij = 4 ν e,i ν e,j (ν e,i + ν e,j ) .
The momentum equation in the variable-h SPH discretization is given by, where, and η = 0.001h 2 i is a small number added to ensure a non-zero denominator in case when i = j.
Remark. We do not employ any artificial viscosity in our benchmark cases. We note that the proposed scheme is not conservative due to shifting, the adaptive-h correction terms, and the non-standard form of the pressure gradient.

Particle shifting
We use a limited form of the particle shifting technique of Lind et al. [32] which is based on evaluating the gradient of the kernel function. A particle with an index i at a current position r i is shifted to a new position r i as, where, where ξ is the point of inflection of the kernel [33], and h ij = (h i + h j )/2. For quintic spline the point of inflection is ξ = 0.759298480738450. We found that using ρ j in the volume approximation makes the shifting less effective and hence have used ρ 0 . We limit the shifting by restricting the movement of particle which is shifted by more than 25% of its smoothing length: We employ shifting while solving the fluid equations and also after our adaptive refinement procedure. Since we use the transport velocity scheme which already accounts for the shifting no additional correction is necessary. However, after the adaptive refinement procedure and subsequent shifting we correct the fluid properties by using a Taylor series approximation. Consider a fluid property ϕ i the corrected value ϕ i is obtained by, The transport velocity is computed using the shifting as,

Boundary conditions
We employ periodic, no-slip, free-slip, no-penetration and the inlet-outlet boundary conditions in our test cases. We enforce periodic boundary conditions by the use of ghost particles onto which the properties are directly copied from the particles exiting the domain through a periodic boundary.
For the no-slip, free-slip and no-penetration boundary conditions we use the dummy particle technique of Adami et al. [34]. Dummy particles placed in uniform layers are used to discretize the wall. The no-penetration is implicitly enforced by using the wall velocity in the EDAC equation [34]. For the no-slip or free-slip we extrapolate the values of velocity of the fluid onto the dummy wall particles by, where the subscript w denotes the dummy wall particles, u i is the prescribed wall velocity, andû is the Shepard extrapolated velocity of the fluid particles indexed by j onto the dummy wall particles i. The pressure on the wall is calculated from the fluid, to accurately impose the pressure gradient, by, where the subscript f denotes the fluid particles, a w is the acceleration of the wall, r wf = |r w − r f |, and h wf = (h w + h f )/2. For the inlet and outlet we use the non-reflecting boundary condition of Lastiwka et al. [35]. First we compute the characteristic properties, referred to as J 1 , J 2 , and J 3 in aforementioned article, of the fluid. Then, we extrapolate the characteristic variables of the fluid onto the inlet and outlet particles using Shepard interpolation. Finally we determine the fluid dynamical properties from the characteristic variables.

Force computation
We compute the forces on the circular cylinder in the flow past a circular cylinder simulation and evaluate the coefficients of lift and drag. Specifically, we compute the forces due to the pressure and the skin-friction on the cylinder by evaluating, which in the variable-h SPH discretization is written as, where the summation index j is over all the fluid particles in the neighborhood of a solid particle indexed by i. We compute the coefficient of pressure drag c d,pressure and skin-friction drag c d,skin-friction , and coefficient of lift c l due to pressure by, where L is the characteristic length of the simulation, U ∞ is the free stream velocity, f visc = j f j,visc and f p = j f j,p is the sum over all the dummy wall particles, and e x and e y are the unit vectors in the x and y directions respectively.

Time integration
We use Predict-Evaluate-Correct (PEC) integrator to integrate the position r i , velocity u i , and pressure p i . The integrator is as follows: We first predict the properties at an intermediate time value n + 1 2 , next we estimate the new accelerations at n+ 1 2 . We then correct the properties to get the corresponding values at the new time n + 1, The time-step is determined by the highest resolution used in the domain, and the minimum of the CFL criterion and the viscous condition is taken:

Adaptive refinement
We first provide a broad overview of the method before delving into the details. The adaptive refinement algorithm involves the following key ideas: • A particle is split if its mass is greater than m max . Note that m max is space varying. The splitting is performed using the approach of Feldman and Bonet [2] and Vacondio et al. [3]. We normally split each particle into 7 child particles in two dimensions.
• A particle i is allowed to merge with another particle j if r ij < (h i +h j )/2 ). The merging algorithm is fully parallel and only particles that are mutually closest to each other are merged. That is, only if particle i's closest allowed merge particle is j and j's closest allowed merge particle is i, will particle i and j be merged. More details on the merging algorithm are provided below.
• When the particles are split they are iteratively merged three times in order to merge any split particles with nearby particles.
• The maximum mass and minimum mass at a particular location are set automatically using a reference m ref parameter that is automatically computed based on the minimum specified resolution and a ratio C r similar to what is done by Yang and Kong [13]. m max = 1.05m ref and m min = 0.5m ref .
• A collection of "background" points is used to adaptively set the minimum and maximum mass values of the fluid particles to control the adaptive resolution. If the body moves, this background is also updated.
The global minimum and maximum size of the particles is specified. Any solid bodies (barring the far-field slip walls) are assumed to be specified at the smallest size. The reference mass increases by the ratio C r from the smallest particle to the largest. This produces a smooth increase in the number of particles in each region.

Adaptive splitting
The algorithm for splitting particles follows that of Feldman and Bonet [2] and Vacondio et al. [5]. If a particle's mass is larger than the maximum allowed mass, m max , then it is split into 7 particles. The original particle is called the parent particle and the split particles are called child particles. Six child particles are placed in a hexagonal arrangement with one child particle at the center as shown in Fig. 1. The parent particle has a smoothing length of h, the child particles have a smoothing radius given by αh and they are placed on a circle of radius h. These parameters α, are normally computed so as to minimize the density error as discussed in [2,5]. We choose the parameters α = 0.9, = 0.4 for the equal mass ratio case. We note that in the present work these are only initial values of the distance and the smoothing length factors. After we split the particles we perform merging followed by shifting and a Taylor series correction. These corrections are accurate due the choice of the values of α and . Subsequently, we use the optimal smoothing length as discussed in the following. The mass of all the particles is the same and is equal to a seventh of the parent's mass. This configuration produces very little error. We use a quintic spline kernel for all computations in this work. The density, velocity, and pressure values of the parent particle are copied to the children. The children also copy the values of the m min , and m max . ǫh h αh Figure 1: Sketch of particle splitting. The parent particle is shown as a dashed circle and has a smoothing length of h. Six child particles are placed in a hexagonal pattern with one child at the center.
We note that in [5], the value of parameter α is 0.9. This implies that the smoothing radius of the child particle is 0.9 times that of the parent despite it having a mass of around a seventh of the parent. Normally in an SPH simulation one tends to choose m = ρ∆x d , where d is the number of spatial dimensions and ∆x is the inter-particle spacing. Furthermore, h = k∆x and k depends on the choice of the kernel. Thus, the value of α = 0.9 is much larger than what one would ordinarily expect. This makes the original approach computationally inefficient and significantly increases the number of neighbors of each particle. This also reduces the accuracy of the method since the smoothing errors are larger. In the present work we find the average mass of particles in the neighborhood of each particle and use this to set the smoothing length using, h = C(m/ρ) 1/d , where C is a constant. In regions where the particle mass is uniform, this attains the ideal h value that would have been set without the use of adaptive resolution. This gives us an optimal h and is therefore computationally efficient. We test the accuracy of our method with a suite of benchmark problems in section 4 and find that this does not affect the accuracy of the method.
The implementation of particle splitting is relatively straightforward. The adaptive splitting may be performed either every iteration or every n adapt > 1 iterations. Any particles whose mass is greater than the m max value are split. Once these particles are identified, the total number of particles that need to be split can be identified. In addition, we also identify the particles that are to be merged as discussed in the next section. Hence, the total number of new particles that need to be created is known. The new child particles are then stored over any unused merged particles and new particles that have been created. Each of these steps are easy to implement in parallel using a combination of elementwise and reduction operations.

Merging particles
The merging algorithm is in principle simple and we use essentially the same approach as discussed in [3]. We note that the smoothing radius of the particles is initially set as discussed in [3]. If we have two particles at locations r a , r b . The location of the merged particle is at, where m m = m a + m b is the mass of the merged particle. The velocity is set using the mass-weighted mean as, A similar form is used for any scalar properties like pressure. The position r m and velocity u m are obtained by ensuring the conservation of momentum. The smoothing radius is set by minimizing the density error [3], where W (x, h) is the kernel function and d is the number of spatial dimensions. Once the entire splitting and merging process is complete, the smoothing length h is set to an optimal value as discussed in the previous section using the average mass of the neighboring particles.
The parallelization of the merge step is however, a non-trivial problem which we discuss here. We wish to use a parallel algorithm that can identify possible merge candidates in one loop over the particle neighbors. The algorithm is designed so each particle can identify a suitable merge partner in parallel. This is achieved using the following approach.
• A particle i is allowed to merge with another particle j if the particle j has not been identified for splitting and ). All neighbors of particle i are searched and the closest particle index closest_idx that satisfies these criterion is identified. This is a completely parallel operation.
• If the i'th particle's closest_idx is j, and if the closest index of the j'th particle is i, then the two particles may be merged. Otherwise the particles are not merged.
• Once a pair of merging particles are identified, the particle with the smaller numerical index value is retained and the particle with the larger index is marked for deletion.
This algorithm is entirely parallel and can be implemented on a CPU or GPU very easily. In fact, these computations may be implemented easily in the context of a SPH calculation. After the identification is complete, one can easily identify the particles that need to be deleted or merged.
The above algorithm may run into pathological particle configurations which will not merge enough particles. However, we find that this does not happen in practice and the algorithm works rather well especially since the particles are constantly moving and are homegenized by the use of a particle shifting procedure. Figure 2: Sketch of particles splitting and then being merged. On the left side is a set of particles that split into 7 children each. On the right, these particles are merged to reduce the number of particles. Note that only one round of merging is complete at this stage.

Merge
It is important to note that when the particles are split, one particle is split into 7 (as discussed in section 3.1). In order to reduce the number of particles we also iteratively perform merging using the same algorithm as discussed above. The reason we choose to split particles into 7 and then merge is that this tends to produce much lower errors since the particle distributions after splitting are more uniform and this makes it more effective to find merge partners. Since the merging is performed pairwise and we desire that on the average each particle be split into two particles, we must have at least three merges. Increasing the number of merges is computationally expensive so we limit it to three. In Fig 2, we show on the left two columns of parent particles that are moving. As they move to the right, they split into 7 children each. These are merged once to produce the particles on the right. With a subsequent merge the remaining small particles are also merged into larger particles depending on the allowed maximum and minimum masses. The figure indicates that the particles are disordered. In order to correct these we perform particle shifting iteratively three times using eq. (18) and correct the properties of the fluid using eq. (21). In subsection 3 we show some particle plots (see Fig. 4) where one can clearly see that the particles are uniformly distributed.

Automatic adaptation
The key part of the adaptive split and merge algorithm is in setting the appropriate m max (r) and m min (r) spatially. In simple cases, it is possible to manually assign the appropriate reference mass for different spatial regions. On the other hand for more complex cases we may not be able to set this manually. For example when simulating the flow past a bluff body, we would like to prescribe the minimum and maximum resolutions and automatically define the reference mass based on the distance from the solid body. In addition when the solid body moves, the reference mass should be suitably updated. Finally, the algorithm should also support solution adaptivity. We first discuss the simpler case of geometry dependent spatial adaptation and then discuss how solution adaptivity can be added.
We setup the discussion in the context of wind-tunnel-like problems where a collection of stationary or moving solid bodies is placed in a stream of fluid with a suitable inlet and outlet. In these class of problems, the solid body typically defines the highest resolution since this is where the largest gradients are observed.
We use the term size of a particle to refer to the inter-particle spacing ∆s. We determine a suitable reference mass in a region, m ref and then set m min = 0.5m ref and m max = 1.05m ref . The size of the particle immediately determines its m ref , for in two-dimensions, m ref = ρ(∆s) d where d is the number of spatial dimensions. In order to smoothly vary the regions, we use a parameter 1.05 ≤ C r ≤ 1.2. The sizes of particles in two adjacent regions are in this ratio, i.e. ∆s k+1 = C r ∆s k , where k indicates a layer of particles with a similar resolution. We assume that the minimum resolution ∆s min and the maximum resolution ∆s max for the particles are known quantities. We note that h = 3∆s for the simulations in this work.
The Lagrangian nature of the SPH method makes it difficult to use the fluid particles to themselves define the reference mass. Instead, we employ a set of stationary background particles. These background particles are not involved in the computation of the governing equations of motion of the fluid or solid. They are merely used to set m ref based on the requirements. The background particles are initially setup with a constant size of ∆s max . The solid geometry of interest is discretized at a resolution of ∆s min . Given these, we initialize the background when the simulation starts as follows.
1. Iterate over all the background particles. If a background particle has a solid particle as a neighbor, then the background particle is marked as being near a boundary. In our implementation, we have a simple integer mask which is set to the value 1. These particles are set to have the smallest size (the same as that of the solid particles). We call these the boundary background particles. 2. Once the boundary background particles are identified, we iterate over the remaining particles and find the minimum (∆s min ), maximum (∆s max ), and geometric mean (∆s avg ) of the sizes of the neighboring particles. If the ∆s max /∆s min < C 3 r , this suggests that the regions are near the ideal distribution, and we set the size of the particle to be equal to C r ∆s min . If ∆s max /∆s min >= C 3 r , then we set the size of the particle to ∆s avg . This allows the particle sizes to be refined rapidly in the initial stages when most of the particles are at the highest resolution. When the distribution is nearing the desired distribution we ensure that the nearby layers are such that ∆s k+1 = C r ∆s k , where k indicates a layer. We note that the size of the particle immediately determines m ref , m min , and m max . 3. Once the reference mass of the particles is set, the particles are split if required using the same splitting algorithm as used for the fluid particles. 4. The smoothing length of the particles is now set such that the number of neighbors is roughly the same. Equation (41) is used for this purpose and is discussed below. The background particles are also moved to distribute them uniformly using the same PST method as used for the fluid. Both the method of [3] or [32] work well although we use that of [32] in this work since it is parameter free and works very well. These two operations of setting the smoothing length and using a PST are repeated three times.
The equation used to iteratively set the smoothing length of the background particles is the same as that used in [13] and is reproduced here, where N r is a reference number of neighbor particles, N n i is the number of neighbors for particle i at iteration level n. This approach ensures that each particle has close to N r neighbors eventually. We reiterate that this algorithm is only used for the background particles so that they smoothly vary and the computations of the references masses are smooth. For the two-dimensional flow problems considered here, we use N r = 48.
In order to initialize the background particles, the above steps are repeated 3 log(∆s max /∆s min )/ log(C r ) times to setup the initial background. The fluid particle resolution is set by finding the minimum of the background particle reference mass in its neighborhood. Thus the background particles only define the spatial resolution for the fluid particles.
In Fig. 3 we show the background particles and the corresponding number of split levels with 0 being the smallest size particles and 7 being the largest. In this case the minimum particle spacing is 0.1 and 0.4 is the maximum spacing. We choose a C r = 1.2 and this generates roughly 8 regions with differing m ref values. The fluid particles created for this distribution of particles is shown in Fig. 4. Here we can see that there are only 4 layers since when particles split they effectively split when the mass from one layer to the next jumps by a factor of two.
When the solid bodies move, the algorithm above is executed once every few iterations. This automatically adapts the reference mass distribution in a smooth fashion. Since the motion of the bodies in each time step is typically quite small and a fraction of the local smoothing length, we only need to perform one iteration of the above.
In Fig. 5 we show the case of two unit square solids placed in a fluid, the background particles are shown and the particle size is smoothly decreasing towards the solid geometry. We move each square by 0.05 units away from each other in each step and update the background by performing the steps discussed above once each time step. We do this 60 times and the resulting background particles are shown in Fig. 6. As can be clearly seen, the background mesh adapts to the moving solid. This shows that the algorithm can comfortably handle moving geometries.
We note that the current method may also be used to setup a specific user-defined region with a desired resolution. This is done by creating a set of particles that serve as a solid body but are only used to set the   boundary background particles. These particles do not participate in any fluid-solid computations. Thus the approach offers a convenient way to define user-specified regions with different resolutions.
While we do not explore this extensively in the current work, it is easy to incorporate solution adaptivity using this framework. Let us assume that there is some solution dependent scalar φ(r) that may be evaluated using the  fluid particles (like the magnitude of the vorticity) and are interpolated onto the background particles. We can use a linear mapping between the range of the values of φ to the minimum and maximum allowed resolution. Once the boundary background particles are identified (step 1 in the algorithm for the background particle) we use the φ value to appropriately set the resolution. The rest of the algorithm then proceeds as before to update the remaining particles.
We show examples of solution adaptivity based on the vorticity, ω in section 4.5. In this case, we compute the absolute magnitude of the vorticity of the fluid particles and interpolate them onto the background particles as the value of φ. Any particles with a value of φ > k max(ω), where k is a user-specified value, are assigned the highest resolution. This approach allows us to track the vorticity adaptively. The approach may be easily extended to use different refinement criterion if so desired. The proposed algorithm can thus handle a variety of different forms of adaptivity.
We want to assess the errors due the splitting and merging. First, we estimate the error in the global density due to merging of particles of different mass ratios at varying separation distance. We consider a square domain Ω of uniformly discretized points and two additional particles which are merging. We compute the error in the global density [5], due to the merging, defined as, (42) where h m is set using eq. (40), subscripts a and b denote particles which are being merged, and subscript m denotes the merged particle. β = ma m b is the mass ratio of the particles being merged. The integration is performed over the region of uniformly discretized points and the particles have full kernel support. Figure 7b shows the error in the global density for two particles merging at varying separation distance and mass ratio. It is clear from the figure that the merge errors are low when the particles are closer and have nearly equal mass.
Next, following Vacondio et al. [5] in order to assess the error due to splitting one particle into 6+1 daughter particles, we compute the error in global density. We consider a square domain Ω of uniformly discretized points and a single SPH particle. We split a single SPH particle into 7 daughter particles, where 6 daughter particles are placed on the vertices of a hexagon centred around the parent particle and one daughter particle at the location of the parent particle. The distance of the 6 daughter particles from the center is controlled by the parameter , where r p = h p , and the smoothing length of the k th daughter particle, h k = αh p , is controlled by the parameter  Figure 7: (a) Error in density by splitting one particle into 7 daughter particles with equal mass. We vary the smoothing length factor α and the distance of the 6 children particles from the parent particle. (b) Error in density while merging two particles of varying mass ratios and separation distance.
α. The global density error is then computed by evaluating, where λ k is the mass ratio of the parent particle to the k th daughter particle, the subscript p denote parent particle, the summation is over all the children particles. The kernel function is computed and the integration is performed over the uniformly discretized points. The parent and all the daughter particles always have full kernel support. We take λ k = 1/7 in this work based on the results of the errors in merging shown above. Figure 7a shows the error in the global density for a particle split into 7 equal mass daughter particles. Based on the results of errors in merging and splitting we choose the values of the daughters' smoothing length factor α = 0.9 and position from the center = 0.4.

Algorithm
In this section we summarize the adaptive resolution algorithm. We start with a given solid body or multiple such bodies that are discretized at the highest desired resolution, with particle spacing, ∆s min . For complex geometries, we may use the particle packing method proposed in [36] to generate uniformly distributed particles for discretizing the solid bodies. We prescribe a coarsest resolution ∆s max as well as the desired C r factor which is typically between the values of 1.05 to 1.2. This effectively determines the width of each refinement layer. One may also manually specify the constraints on the mass in different spatial regions. Finally, we are given a domain of interest; in the problems considered in this work, the domain size is fixed and known a priori. Given this, the algorithm proceeds as follows.
1. We first initialize the background particles as discussed in section 3. 2. Using the background particles, we initialize the fluid particles at the initial time. This is done by splitting and merging the particles based on the reference mass of the background particles. This is discussed in detail in sections 3.1 and 3.2. The particle shifting algorithm of Lind et al. [32] is applied at each stage to get a smooth distribution of particles. 3. The initialized fluid particles along with the given solid particles are then used to simulate the governing equations using an appropriate scheme. In the present work we use a highly modified EDAC-SPH scheme as discussed in section 2. 4. At the end of every iteration, we find the nearest fluid particles to the background particles and set the reference mass of the fluid particles. 5. The fluid particles are adaptively split and merged every n adapt iterations as discussed in sections 3.1 and 3.2. Typically we choose n adapt to be between 1-10. Similarly, the background particles are updated (to accommodate moving bodies or solution adaptivity) every n bg iterations. When the adaptation is entirely spatial and the solid boundaries do not move, the background does not need to be updated at all; we usually set this to 100 or 500 iterations. This parameter is adjustable depending on the requirements. 6. After merging, we shift the particles three times using eq. (18) and correct the properties of the fluid eq. (21). 7. The smoothing length of all the particles is set using the average mass of the neighbors as discussed in section 3.1.
After the initialization is complete, the time-marching procedure of the SPH scheme is started. Before every iteration, the following is performed: (i) we update the background particles every n bg iterations; (ii) thereafter, every n adapt iterations, we split, merge, shift and correct the properties of the particles; (iii) we update the smoothing length of the particles based on average mass and update the nearest neighbor search algorithm to fetch new neighbors; (iv) finally, we execute each time-step of the time-marching scheme as discussed in the section 2.4.
In order to assess the accuracy of the algorithm, we consider the following test case. We consider a square region of side two units on which we discretize the function f (x, y) = sin(2πx) + cos(2πx) using particles of equal size. The inner square region of unit side is split into particles. We consider two cases, the first using the formulation of Vacondio et al. [5] where each particle in the inner region is split into 7 daughter particles, and the second using the proposed method. We then compute the error in function and gradient approximation figs. 8 and 9. The results are summarized in table 1. As can be seen the proposed algorithm is accurate and has far fewer neighbors than Vacondio et al. [5].

Results and discussion
We apply the adaptive resolution technique proposed in this work to the test cases shown below. We first apply our method to the classical numerical test cases with varying Reynolds numbers and compare with established results in the literature. We then simulate the flow past a circular cylinder at Reynolds numbers 40, 550, 1000, 3000, and 9500. We show the details that typically require a large number of particles to capture accurately. We compare the results of our method to the high resolution vortex method results of Koumoutsakos and Leonard [37], and Ramachandran [38]. We also show the results of solution-based dynamic particle resolution, where the particles are adaptively resolved to the highest resolution based on the magnitude of vorticity in the flow, for the flow past a circular cylinder and the flow past a C-shape at Re = 2000. Every figure presented in this manuscript is automatically generated by an automation framework [24]. The open-source code is available at https://gitlab.com/pypr/adaptive_sph.

Taylor-Green vortex
The Taylor-Green problem is a widely used benchmark to study accuracy in SPH [5,8]. The exact solution for the Taylor-Green problem is given by, v + = U ∞ e bt sin(2πx + /L) cos(2πy + /L), where the superscript (.) + indicate dimensional quantity, U ∞ = 1 m/s, b = −8π 2 /Re, Re = U L/ν, and L = 1 m. In the results we use the dimensionless velocities u = u + /U ∞ and v = v + /U ∞ , pressure p = p + /ρ 0 U 2 ∞ , and distance  We consider two Reynolds numbers 200 and 1000, and simulate the problem using the adaptive algorithm proposed in this paper with three different minimum resolutions L/∆x max of 50, 100 and 150. We compare the results with the exact solution and the non-adaptive simulation with a resolution that matches the minimum resolution of the adaptive case. For the Re = 200 case we also compare with the Adaptive Particle Refinement (APR) results of [8]. Figure 10 shows the velocity magnitude and the pressure particle plots for Re = 200 at t = 2.5 s. The velocity and pressure contours show less decay than [8]. Specifically, note that the high velocity regions near the interface of the layers are well maintained in the present work. We show the results for Re = 1000 in fig. 11. It can be seen that the results are as expected and the contours do not show any decay or change in shape. We note that there is a drift in the average pressure and subtract it from the pressure contour plots. This drift in pressure is due to errors in the total volume conservation [27], which is computed as follows, Figure 12 shows the error in total volume. We compare the adaptive and nonadaptive cases with the non-adaptive δ + -SPH results by [27] who simulated the problem with Re = 1000 at a resolution L/∆x = 400. As can be seen in section 4.1, the errors in the present work for L/∆x = 150 are almost five times smaller than the case with L/∆x = 400 of [27]. However, the adaptive cases have a larger errors which drop as the resolution increases.  Figure 13a shows the spatial distribution of the smoothing length h. As can be seen the smoothing length is almost constant in the interior of the respective regions. At the interface between the two regions having different mass the value is changing gradually. Figure 13b show the distribution of number of neighbors of each particle. It can be seen that in the interior of the regions the value is around 30. Whereas, in the interface between the two regions it is larger as would be expected. Since bulk of the particles have minimum number of neighbors the method is efficient.
In fig. 14 and fig. 15 we plot the maximum velocity decay and the L 1   error in the velocity for Re = 200, and Re = 1000 respectively with different minimum resolutions. The maximum velocity decay shows good agreement with the exact solution. We also compare with the non-adaptive case at different resolutions. Although we do not expect greater accuracy than the non-adaptive case due to the presence of lower resolution regions, we expect the errors to be of the same order as that of the non-adaptive case. The L 1 norms reveal that the errors in the adaptive case are almost 2 times the errors in the non-adaptive case for Re = 200, whereas for the Re = 1000 case, the errors are lower. The increase in the L 1 error for the adaptive case at Re = 200 is due to the effect of having an approximately constant number of neighbors in the adaptive region. Normally, the discretization error due to the viscous operator reduces when the number of neighbors is increased. In addition, this error is significant at low Reynolds numbers and hence causes an increase in the L 1 error. This is consistent with the findings of [39] (figure  11 and table 5), and [40] ( figure 15). A theoretical proof of this is available in the work of [41]. At Re = 200 the viscous error dominates, however, at Re = 1000 the viscous effect is not dominant and the L 1 error behaves as expected; the adaptive method has a smaller L 1 error than the non-adaptive case.
In fig. 16 we show the kinetic energy decay for Re = 200 and Re = 1000 at different minimum resolutions. We compare with the exact, non-adaptive, and the APR simulation of [8]. Our results match the non-adaptive decay and barring a slight increase the exact decay, whereas the APR scheme shows comparatively large decay. Figure 17 shows the error in the kinetic energy versus time. For the Re = 200 case the results show an anomalous behavior as seen from the L 1 error. The Re = 1000 cases behave as expected.    The above results show that the proposed method is accurate, displays less dissipation than other recent techniques proposed for adaptive resolution, and requires minimum number of neighbors for bulk of the particles. This makes the proposed method both accurate and efficient.

Gresho-Chan vortex
Gresho-Chan vortex [42] is a two-dimensional inviscid numerical test case with periodic boundary conditions in both the x and y directions. This problem tests, to name a few, the numerical stability of the method, and the conservation properties. Considered as a difficult test case [43], this test case is widely used by the astrophysical community [44,45,46]. The problem is of a rotating vortex inside a domain of unit length where the centrifugal force due to azimuthal velocity balances the pressure gradient. The distances are non-dimensionalized by the domain length L. The initial radial velocity is zero, and the azimuthal velocity in non-dimensional form is given by, where R = 0.2L, and r is the distance from the centre of the vortex located at the origin (0, 0). The artificial speed of sound is c s = 10 m/s, and the smoothing length factor, h/∆x = 1.0. The non-dimensional pressure, balanced by the centrifugal velocity, is given by, where p 0 = 5 is the reference pressure. We adaptively refine a semi-circular region of radius 0.45 around the origin with particles of mass around 0.5 times that of the outer particles. We simulate the problem for t = 3 s. We compare our results with the exact solution and the non-adaptive cases. We consider two different minimum resolutions L/∆x max of 50 and 100. The particles are initially placed on a uniform Cartesian grid. Figure 18 shows the particle positions at t = 3 s. It is difficult to assess the difference between the simulations from this result. Figure 19 shows the magnitude of the velocity of all the particles in the domain as a function of the distance r from the centre of the vortex. The red-line indicates the exact velocity magnitude. The figure also indicates the L 1 norm of the error, which is computed as, The plot shows decay and noise in the velocity magnitude. This is an inviscid problem and our simulation does not employ any artificial viscosity. It is therefore highly sensitive to small perturbations. The results show that the particle splitting and merging process introduce a small amount of noise in the simulation. However, the results show that this is only slightly dissipative. One can also see that these results are as good if not better than the existing results [43,46]. Further, Hopkins [46] mentions that splitting and merging can be noisy and diffusive. However, the present results show that careful splitting and merging of particles as we have done produces acceptable results. Figure 20a shows the angular momentum of the system as a function of time. For the non-adaptive cases we can clearly see a small amount of dissipation which reduces as we increase the resolution. A similar trend follows for the adaptive cases. Figure 21 shows the evolution of linear momentum in x and y directions where higher resolution has better conservation. But due to the reasons mentioned in the remark of section 2 the conservation of momentum does not hold. We note that the mass is exactly conserved in all our cases. The adaptive L/∆x max = 50 case is more dissipative compared to the L/∆x max = 100 case. To further study this, the maximum velocity evolution over time is shown in fig. 20b. In this figure we see that for the adaptive L/∆x max = 100 case, the maximum velocity does not decay significantly. These results affirm the accuracy of our adaptive algorithm.

Two dimensional lid-driven cavity
Lid-driven cavity is a two-dimensional viscous problem with solid boundaries. We study this problem with two different Reynolds numbers 100 and 1000. We compare our results to those of Ghia et al. [47]. The domain length L is 1 m, and the top wall is moving with a velocity U wall of 1 m/s. We consider a square domain  fig. 22. The artificial speed of sound is c s = 10 m/s, the smoothing length factor, h/∆x = 1.0, and the reference density ρ 0 is 1 kg/m 3 . We non-dimensionalize the velocity by the wall velocity, u = u + /U wall and v = v + /U wall , the pressure as p = p + /ρ 0 U 2 wall , and the lengths by the domain length L.  Figure 22 shows the velocity magnitude distribution and pressure distribution for Re = 1000. We use 3 layers to simulate this problem. The outer layer of particles are at the highest resolution with the particle mass corresponding to a resolution of L/∆x min = 150. The middle region is at twice the mass of the outer region, this corresponds to a resolution of L/∆x = 100. The inner most resolution is the coarsest of all with an effective resolution of L/∆x max = 75. Figure 23 shows the centerline velocity profiles at Re = 100. The results match well with the results of [47]. In fig. 24 we show the centerline velocity profiles for Re = 1000. It can be observed that as the resolution is increased the centerline profiles show a good agreement with [47].

Flow past a circular cylinder
We study the flow past a circular cylinder problem for five different Reynolds numbers ranging from 40 to 9500. We plot the coefficients of pressure drag and skin friction as a function of time. We compare the results with the high resolution vortex method of Koumoutsakos and Leonard [37], and Ramachandran [38]. Figure 25 shows the domain setup, where the diameter of the cylinder D = 2R = 2 m. We use a non-dimensional time T = tU ∞ /R. We initialize the flow at T = 0 with the potential flow solution.
The inlet velocity is 1 m/s, and the solid walls are inviscid. In order to minimize the reflection of the initial, undesirable, pressure waves from the walls we employ the non-reflection boundary conditions of Lastiwka et al. [35] on the inviscid walls.
We simulate the problem with fixed refinement zones up to T = 6. For all the simulations, the coarsest resolution in the domain is D/∆x max = 4. We vary the finest resolution D/∆x min from 100 to 500. For the Re = 9500 case we use a finest resolution D/∆x min = 1000. Unless explicitly mentioned we use a C r value of 1.08 for all the problems. The parameters used in these simulations are summarized in   First, we demonstrate the advantages of using adaptive particle refinement over the non-adaptive case. We use two Reynolds numbers 1000 and 3000.
We simulate the problem up to T = 6 with both the non-adaptive case, using a resolution D/∆x = 50, and the adaptive case, using two different maximum resolutions D/∆x min = 50 and 100. We use solution adaptivity based on the vorticity and this aspect is explored in greater detail in the next section. Figure 26 shows the coefficients of pressure drag for the Re = 1000 case. The coefficients match closely for both the non-adaptive and the adaptive cases, the adaptive case with D/∆x min = 100 is slightly better in comparison. The differences between the adaptive and non-adaptive are not easy to assess in this case. On the other hand, for Re = 3000 the advantage of using the adaptive resolution is clearly seen in fig. 27 and summarized in table 4. As can be seen, the adaptive and non-adaptive cases match at D/∆x min = 50. Although, at D/∆x min = 50 the adaptive case uses 35 times fewer particles and is 22 times faster than the non-adaptive case. Given the efficiency of the adaptive particle refinement, we increase the resolution in the adaptive case to D/∆x min = 100 and further to 200 and observe significantly better results in fig. 27. The results match those of the vortex method quite closely. For the D/∆x min = 200 case the number of particles is 2.6 times that of the D/∆x min = 50 case with the adaptive particle refinement. This simulation still requires 13 times fewer particles than the non-adaptive case at a much lower resolution of D/∆x min of 50. These results clearly indicate the importance and performance of the adaptive particle resolution.    numbers. We only show the results where the finest resolution is 500. The results are in good agreement with that of [37,38]. For the case of Re = 9500 the results of [37] predicts slightly (note the logarithmic scale) higher skinfriction drag, but our results match closely to that of [38]. In figs. 29 and 30 we plot the coefficient of pressure drag for the Reynolds numbers, 40, 500, and 1000, 3000 respectively. Our results differ from the established results at the start for up to T = 0.5. This is due to the weaklycompressible nature of our flow for where an initial pressure wave is required to set the velocity from the potential start to the viscous profile, whereas the established results use incompressible flow. Thereafter, our results match closely with increase in the maximum-resolution.  For the Re = 9500 case shown in fig. 31 we use the finest resolution D/∆x min = 1000 with a Cr = 1.15. The change in C r is due to time and computational constraints. This is the highest resolution used in our simulations. Even though the characteristic features of the drag coefficient profile match with the established results the curve does not reach the maximum, the trends are consistent with the established results. We note that the resolution used by [38] corresponds to a finest resolution of D/∆x min = 1250, [37] use a million vortices for their simulation, and the present simulations employ around 200,000 fluid particles in the entire domain. We would also like to note that after three seconds maintaining symmetry is difficult and even with the DVH results of Rossi et al. [48] there are significant differences in the drag force. The present results are clearly in good agreement given the variation in literature.  Figure 32 show the radial velocity along the axis of symmetry on the rear side of the cylinder for the Reynolds numbers 3000, with D/∆x min = 500, and 9500, with D/∆x min = 1000 and Cr = 1.15, at different times. In both the Reynolds numbers the results are in good agreement with Ramachandran [38] while the results of Subramaniam [49] show slight difference for larger time in the Re = 3000 case. The results are compared with Ramachandran [38] and Subramaniam [49]. Figure 33 compares the proposed formulation vorticity distribution with Durante et al. [50]. The simulation is run at the Reynolds number 1000, and times T = 6.4, and 12.8 are shown. In the adaptive SPH figure (left) the particles are sized proportional to the mass. There are some differences in the color as a slightly different color map was used. The distribution show a good similarity.  [50]. Adaptive-SPH particles' size is proportional to their mass. The finest resolution, D/∆x min , is 200, and the coarsest resolution, D/∆x max , is 40. The value of C r is 1.08.
In fig. 34 we compare the vorticity distribution at the Reynolds number 9500 with Ramachandran [38]; times T = 1, 2, and 3 are shown. There are some differences in the color as a slightly different color map was used in [38]. The vortices appear to maintain the symmetry, and the secondary, tertiary, and further vortices generated at the boundary layer are captured well. The boundary layer at the leading edge of the cylinder is clearly observed. As the simulation progresses, the vortices grow big and move across different layers having different smoothing lengths since there is no solution adaptivity used in this case. At T = 3, the plot clearly shows the primary vortex at a different resolution than the boundary layer.   Figure 35 shows a histogram of the number of neighbors in the overall simulation at T = 6 for the resolution D/∆x min = 1000. This shows that for a majority of particles the number of neighbors are at 30. The highest number of neighbors in the simulation is at 42. This shows the optimal neighbor distribution further maximizing the performance. In fig. 36 we show the smoothing length distribution for the same case. The left side shows the whole domain and the right is a zoom-in near the cylinder. The smoothing length varies across a large number of scales by a factor of 250. Even near the cylinder it varies by about a factor of 20.

Solution adaptivity
In this section we demonstrate the solution-based adaptivity. We consider a flow past C-shaped body at Re = 2000 and compare our results with Rossi et al. [48], and Sun et al. [11]. The domain dimensions are given in fig. 37, and the smoothing length factor, h/∆x = 1.2. The outer diameter D is 1 m. We first perform the simulation without adaptivity and compare our results. The minimum-resolution D/∆x min is 25 and the maximum-resolution D/∆x min is 200. The minimum and maximum-resolution match the respective resolutions of Sun et al. [11]. We simulate the problem for T = 30. In fig. 38 we show the coefficients of total drag and lift. The results are in good agreement with Rossi et al. [48]. The initial noise within T = 1 is due to the weakly-compressible nature of our formulation. We compare the number of particles used in our simulation with the simulation of Sun et al. [11]. Sun et al. [11] does not mention the total length of the domain instead provides the dimensions [2.75, 5.75] × [−0.8, 1.2] of an inner rectangular domain containing the C-shape body, where the minimum-resolution is 50. We estimate the number of particles inside this domain to be approximately 98,000, whereas for our simulation the number of particles in this domain is about 38,735. This shows that we use 2.53 times lower number of particles and achieve significantly better results. This further demonstrates the efficiency and the accuracy of our method. Now, we simulate the flow past C-shape using solution-adaptivity, where the vorticity in the flow is monitored and particles with absolute vorticity value above 5% of the maximum vorticity are resolved to the highest resolution. In this simulation we use the maximum-resolution D/∆x min of 125. The simulation is performed for T = 20. In fig. 39 the coefficients of total drag and lift are shown which are in good agreement with the results of [48]. Figure 40 shows the vorticity distribution of the particles for the nonsolution adaptive and the solution adaptive cases at T = 20. It can be seen that for the non-solution adaptivity case the trailing vortices are not refined after they move certain distance away from the C-shape body, whereas in the solution-based adaptivity the trailing vortices are resolved to the highest resolution, based on the cut-off criteria stated above. In fig. 41a we show the zoomed-in view near the C-shape body for the solution-based adaptive case,  [48] and Sun et al. [11]. and in fig. 41b we show the smoothing length distribution demonstrating the effect of solution adaptivity. The results of this section demonstrate the accuracy and efficiency of the proposed method even when there is a large change in the resolution. The method shows little dissipation and it is capable of performing a highresolution simulation to capture all the features of a Re = 9500 flow. The results show good accuracy with the solution-based adaptivity.

Conclusions
In this work we have proposed an accurate and efficient method to handle adaptive resolution in the context of weakly-compressible SPH. This is achieved using (i) an accurate EDAC scheme [21] along with the recent corrections of [26], the use of variable-h corrections of [5], and particle shifting [32]; (ii) adaptive splitting and merging of particles where care is taken to ensure that the number of particles is minimum and the number of neighbors is optimal. We employ background particles to specify the regions of refinement. Importantly, the method allows for specifying fixed regions of refinement, automatic geometry-based refinement, and automatic solution-based adaptivity elegantly in the same framework. The algorithms employed are parallel. We provide an open-source implementation of the entire algorithm along with complete automation of all the results presented in this work.
We demonstrate the accuracy of the method using several benchmarks. The Taylor-Green and Gresho-Chan benchmark problems clearly demonstrate that the method is not diffusive and is more accurate than other recent adaptive refinement techniques.
We perform simulations at unprecedented resolution for the flow past a circular cylinder for a variety of Reynolds numbers in the range 40 to 9500. For example, at Re = 9500 we use a resolution of D/∆x min = 1000 and D/∆x max = 4 giving a ratio of length scales of 250. This requires 16 levels of refinement with a domain size of 25D × 15D, requiring only 200,000 particles. The results are in good comparison with that of [37,38] who also employ similar number of particles. This shows the effectiveness and accuracy of the adaptive resolution method.
For a Reynolds number of 3000 with a D/∆x min = 50 we are able to obtain similar accuracy with the adaptive refinement using 30 times less particles, with a 25-fold speed improvement when compared with that using a fixed resolution.
The current work has focused on the weakly-compressible SPH method. We have demonstrated the method for two-dimensional problems without a free-surface. The method in principle should work with a few modifications for free-surface problems, as well as three-dimensional problems. We plan to explore the adaptive particle refinement applied to incompressible SPH and three-dimensional problems in the future.