6.2.1 Finding the redistribution amounts

The key to the redistribution method is to find a positive solution to the system of linear equations (4.8) and following for the redistribution fractions $f^n_{ij}$. The system is linear, but usually not square: the number of unknown fractions $f^n_{ij}$ is the number of vortices in the neighborhood, while the number of equations is determined by the order of accuracy $M$ desired. In this subsection we will discuss our strategy for obtaining a positive solution for the $f^n_{ij}$, assuming that one exists. The question what to do if no positive solution exists will be addressed in subsection 6.2.3.

The problem of finding a nonnegative solution to an underdetermined system of equations is the standard `phase I' problem in linear programming that can be solved by slack variables. However, following Van Dommelen [235] we will use a different approach. First, we note that the fractions $f^n_{ij}$ must be in the range [0,1]. We may shift the origin to the center of that range, by defining

$$w_j \equiv f^n_{ij} - {\textstyle\frac12} \ ,$$ (6.4)

where the additional dependence of $w_j$ on the vortex $i$ and the time-step $n$ is to be understood. In terms of the $w_j$, a solution is acceptable if the maximum norm of the solution vector, $\Vert\vec w\Vert _\infty \equiv \max_j\{\vert w_j\vert\}$ is less than or equal to $\frac12$.

Our approach is to find the solution for $\vec w$ with the least maximum norm. If the maximum norm is less than or equal to $\frac12$, an acceptable solution has been found. On the other hand, if the maximum norm exceeds $\frac12$, it must mean that no acceptable solution exists. In that case we create more vortices as described in subsection 6.2.3.

The least maximum solution algorithm used in our computations is described in the next subsection. We did do some comparative testing of this algorithm against a standard library routine (IMSL) for the phase I linear programming problem. We found that the number of iterations in the methods was about equal, but that the library routine ran about two times more slowly, possibly due to the extensive safeguards in its implementation. It appears that computational speed is not an important consideration in selecting the method.

However, the least maximum procedure will create a strictly positive solution if one exists, while the linear programming method for the phase I problem will select the minimum number of vortices for the redistribution. As a result, the least maximum procedure tends to spread out the vorticity somewhat better.

In this study, the least maximum problem is solved from scratch for each vortex at each time-step (even for the Stokes flow in which a single solution could have been used for all time-steps). This is a very inefficient approach, since the systems are almost unchanged from one time-step to the next. The relative locations of the vortices in the neighborhood change only by an amount $O(h_v^3)$ during a time-step. This means that a single solution can be used over an asymptotically large number of time-steps. Additionally, in our time splitting we perform two diffusion steps back to back. We solve each from scratch although they are identical.

The disadvantage of using a single solution over many time steps is that some information has to be stored from one time-step to the next. For example, the fractions $f^n_{ij}$ could be stored and updated at each time-step until one turns negative, at which time the system could again be solved from scratch. Alternately, we could merely store the information what fractions have magnitude less than the maximum norm. This is sufficient information to solve a least maximum problem quickly. In any case, our computational times for the redistribution method can presumably still be improved significantly.

We did use one shortcut in our procedure. As a preconditioning to finding the least maximum solution, we performed a Gram-Schmidt orthogonalization on the rows of the system. This orthogonalization directly determines the least length solution, and we found that in about 60% of the cases, the least length solution was positive. Thus we could skip the determination of the least maximum solution in the majority of cases.

There are also tests that could be performed to decide a priori that a system has no acceptable solution: according to estimates given in Appendix A, there must be at least one neighborhood vortex at a distance of more than $\sqrt{4\nu\Delta t}$, the maximum horizontal and vertical distances should be at least $\sqrt{2\nu\Delta t}$, and at least $4\sqrt{\nu\Delta t}/(R+\sqrt{R^2+8})$ in any direction. For third-order accuracy or higher, there should be at least one vortex within a distance $\sqrt{8\nu\Delta t}$.