A fast, matrix-based method to perform omnidirectional pressure integration

In order to convert the explicit omnidirectional scheme into an implicit formulation, we first examine how the line integrals are constructed. For a discrete integration, the pressure at each new point, P_C , is computed from the pressure at the adjacent point, P_a , plus the pressure gradient along the path, $\vec r$.

$$\begin{align} P_C = P_a + \int ^ {\vec x_C} _{\vec x_a} \frac{\partial P}{\partial \vec r} \cdot \textrm{d}\vec r. \end{align} \tag{ 3 }$$

The collection of discrete cells with pressure gradient information (computed from velocity field data) that make up a path across the domain can be defined in multiple ways (figure 1). In this work, we will consider two approaches, shown here in 2D for simplicity. Both approaches considered assume a rectangular, evenly-spaced grid. Irregular grids will not be considered throughout this work.

Zoom In Zoom Out Reset image size

The first type of path, herein dubbed 'Face-Crossing scheme' (figure 1(a)), is made up of all cells that a given path $\vec r$ crosses, and all integrations are made across connected cell faces. The integral path starts from a boundary point B and extends to one of the interior cells C. The second type of path examined, dubbed 'Cell-Centered scheme' (figure 1(b)), consists of the cells with the closest centers. In this case, integrations may be across adjacent cell faces or diagonals. In 2D, if the ray angle θ is steeper than the cell diagonal (i.e. $\theta \gt \textrm{atan}(\Delta y/\Delta x)$), then the cell-centered integration path passes through each y coordinate once and the x coordinate is rounded to the nearest cell for a given ray. Conversely, if the ray angle is less steep than the diagonal ($\theta \lt \textrm{atan}(\Delta y/\Delta x)$), then the path passes through each x coordinate once and the y coordinate is rounded instead.

For the full omnidirectional pressure calculation, the pressure in any cell C will be found from the average of all line integrals that cross the cell (figure 1(c)), which we will determine using the parallel-ray omnidirectional approach (Liu et al 2016). In that method, line integrals of pressure are found for closely spaced rays separated by a shift s and oriented along multiple angles θ.

The proposed method extends the rotating parallel ray approach presented by Liu et al (2016) by taking the limit where the number of line integrals is infinite. Thus, we consider all possible ray displacements s and all possible ray angles θ. As an example, let us analyze the specific case where the rays coming from the east cell E are being counted, as depicted in figure 2.

Zoom In Zoom Out Reset image size

If the cell E is an interior cell (i.e. none of the rays at any angle θ come from a boundary), then the rays can be 'counted' through the integration of a connectivity variable A_E :

$$\begin{align} A_E = \int_{-\pi}^{\pi}\int_{-\infty}^{\infty}R_{CE}\left(s,\theta\right) \textrm{d}s \textrm{d}\theta \end{align} \tag{ 4 }$$

where $R_{CE} = 1$ if C and E are connected by a path given some s and θ; and 0 otherwise. The inner integral in s can be replaced by the projection of the face connecting C and E (in 2D, the interface is an edge), which has height Δy in the direction θ. Also for the E cell, the valid integration bounds for θ are $-\pi/2 \lt \theta \lt\pi/2$, and the integral A_E becomes:

$$\begin{align} A_E^{n+1} = \int_{-\pi/2}^{\pi/2}\Delta y \text{cos} \theta \textrm{d}\theta = 2\Delta y \; \; \textrm{(Interior Cell)}. \end{align} \tag{ 5 }$$

Note in figure 2(a) that all rays coming from the cell E implement the same equation: $P_C^{n+1} = P_E^{n+1} - (\Delta x/2) [\,f_x(C) + f_x(E)]$, where $f_x(j)$ is the x component of the source function evaluated at the cell j and the pressure $P_j^{n+1}$ represents the pressure of the cell j evaluated at the $(n+1)$th iteration of the algorithm. The discretization of the source term $\vec{f}$ in this implementation is a second-order accurate approximation at the cell face, which is an averaging between the values at the center of the cells C and E. For an interior cell, equation (5) defines the ray count $A_E^{n+1}$, whereas $A_E^{n} = 0$ because no pressure values from the previous iteration are required.

However, when the cell E shares one or more faces with a boundary, then the rays intersecting the boundary faces will effectively begin at the cell E, meaning cell E is the boundary cell where the integration starts. Therefore, the pressure value used for cell E has to come from the previous (nth) iteration, as depicted in figure 2(b). For the 'zeroth' iteration, the values can be initialized to zero (Liu and Moreto 2020) or any other arbitrary value. Since for the boundary cell case some rays use pressure values from the $(n+1)$th iteration and some rays use values from the nth iteration, the ray count $A_E^i$ (where i is the iteration number) will be different than equation (5). Depending on the location of cell C, the values of $A_E^i$ will vary. As an example, the values for the boundary topology depicted in figure 2(b) are:

$$\begin{align} \begin{aligned} A_E^{n} = \int_{0}^{\textrm{atan}\left(\Delta y/\Delta x\right)\Delta x} \text{sin} \theta \textrm{d}\theta + \int_{\textrm{atan}\left(\Delta y/\Delta x\right)}^{\pi/2} \Delta y \text{cos} \theta \textrm{d}\theta \\ A_E^{n} = \Delta x + \Delta y - \sqrt{\Delta x^2 + \Delta y^2} \end{aligned} \end{align} \tag{ 6 }$$

$$\begin{align} A_E^{n+1} = 2\Delta y - A_E^{n} = \Delta y - \Delta x + \sqrt{\Delta x^2 + \Delta y^2}. \end{align} \tag{ 7 }$$

Given the forms of the integrals, the 'ray count' $A_j^i$ is effectively the average of edge length (Δx or Δy) projected at all angles θ that come from cell j and use the pressure value from iteration i. For the 'face-crossing' scheme, we can repeat the application of equation (5) to the other faces to find the total projection area from all rays at all angles to be:

$$\begin{align} A_{tot} = \sum_j \sum_i A_j^i = 4\left(\Delta x + \Delta y\right) \end{align} \tag{ 8 }$$

since according to equation (5), $(A_E^{n}+A_E^{n+1}) = (A_W^{\,n}+A_W^{\,n+1}) = 2\Delta y$ and similarly $(A_N^{n}+A_N^{n+1}) = (A_S^{n}+A_S^{n+1}) = 2\Delta x$ due to the symmetry of the problem.

We can then define an implicit equation that implements the averaging process from all rays by performing a weighted average of the finite-difference approximations of equation (3) from all adjacent cells:

$$\begin{align} P_C^{n+1} = \sum_i \sum_j w_j^i P_j^i - \sum_{j} \Delta j \left(w_{j}^{n+1} + w_{j}^{n}\right) \frac{f_j\left(j\right)+f_j\left(C\right)}{2} \end{align} \tag{ 9 }$$

where $P_C^{n+1}$ is the updated value of the pressure at a given cell C at the next $(n+1)$th iteration, and the iteration i can take the values $i = \{n,n+1\}$. The cell j can take the values $j = \{C,E,W,N,S,F,B\}$ for the 'face-crossing' scheme, representing the cardinal directions East, West, North, South, Front and Back, respectively. For the cell-centered scheme, the cardinal directions can also take diagonal values $j = [\{p\},\{pq\},\{pqr\}]$ representing all 27 cells in a 3×3 cube surrounding cell C, including itself. j is also a directional subscript, so $\Delta E = \Delta x$ and $\Delta W = -\Delta x$, and similarly for the other orthogonal and diagonal directions. The source function f_j (at direction j) in the last summation is averaged between cells j and C to approximate the source value at the corresponding face, following a second-order accurate approximation.

Equation (9) simply describes the ray averaging process in the omnidirectional integration approach. The weights $w_j^i$ represent the fraction of the rays that implement each sub-equation (equation (3)), utilizing the values of pressure from the appropriate neighboring cells depending on which direction the ray is coming from. These weights are defined as:

$$\begin{align} w_j^i = \frac{A_j^i}{A_{tot}} \end{align} \tag{ 10 }$$

except for $w_C^{n+1}$, which is defined as:

$$\begin{align} w_C^{n+1} = 0. \end{align} \tag{ 11 }$$

This is because $P_C^{n+1}$ is already counted in the left-hand side of equation (9). The remaining weights are computed from the ray integrals $A_j^i$ as previously discussed.

Note that equation (9) can be defined for all cells in the domain, forming a sparse matrix form:

$$\begin{align} \left[W^{\,n+1}\right] \left\{P^{n+1}\right\} = \left[W^{\,n}\right] \left\{P^{n}\right\} +\left\{S\right\} \end{align} \tag{ 12 }$$

where $[W^{\,i}]$ are weight matrices and $\{P^{i}\}$ are the pressure vectors at iteration i. The source term vector $\{S\}$ represents the last summation in equation (9) and is computed only once based on the velocity measurements from PIV. For a minimal example of how equation (12) would be constructed for a 3×3 domain, the reader is directed to the appendix of this manuscript. Furthermore, the weight matrices $[W^{\,i}]$ can be precomputed for a given grid topology and can include the effect of cells with missing data (say, from regions not illuminated by the laser or containing occlusions/reflections). A method to consider any topology of missing cells will be detailed in the following sections.

The solution method proposed follows the outline provided in algorithm 1. This algorithm was implemented in Matlab as a CPU code (only 2D) and in CUDA-C++ as a GPU code (2D and 3D). For the Matlab implementation, the weight matrices $[W^{\,i}]$ are stored as sparse matrices and equation (12) is solved using Matlab's default linear equation system solver (mldivide or 'backslash' operator). In the GPU CUDA-C++ implementation the Conjugate Gradient (CG) method is used to solve equation (12).

Algorithm 1. Outline of the matrix-based omnidirectional integration algorithm.
Require: Source field $\vec{f}$ from PIV data
Compute: Source vector $\{S\}$ (equation (9), terms using $\vec{f}$)
Compute: cjk and $w_i^{\,j}$ (details in section 2.2)
Compute: Weight matrices $[W^{\,i}]$ (equation (10))
Initialize: $\{P^{n = 0}\} = 0$;
while Residual $\gt$ Threshold do
Solve: equation (12) for $\{P^{n+1}\}$
Compute: Residual (equation (40) or (41))
$\{P^{n}\}\gets \{P^{n+1}\}$
end while

For the GPU implementation, a total of 16 scalar fields with the mesh size are required (7 for the CG solver, 6 for the weight matrices, 3 for the source field $x,y,z$ components). Stored as doubles, the code requires $16\cdot 8N$ bytes of GPU memory, where $N = N_x N_y N_z$ is the total number of cells in the domain. On the Nvidia RTX A4000 GPU used in this work, which has 16GB of memory and retails at ∼$1000 at the time of writing, a mesh size of at most $N = 512^3 = 134$ million cells can be computed. Since the CG solver could potentially compute the weight matrices on-the-fly, it is possible to trade off computational time to increase the size of the domain, requiring only 10 scalar fields to be stored in memory. This possibility was not explored, since a domain of 134 million cells was deemed sufficient for most Tomographic PIV applications currently being considered by the experimental fluids community.

When working with experimental data from PIV/PTV, it is common to have cells in the 2D or 3D grid that are masked out due to missing data. Missing PIV vectors can occur when laser reflections from a model surface are present in the field of view; when a model surface is occluding the view of the flow; or when the laser sheet is not sufficiently expanded to illuminate the entire field of view of the cameras. Alternatively, in Stereoscopic PIV and Tomographic PIV there can also be missing vectors where the views from all cameras do not overlap. A typical flow field with vectors masked out due to the presence of an aerodynamic model is exemplified in figure 3(a).

Zoom In Zoom Out Reset image size

Vectors could also be missing inside the field of view due to post-processing or thresholding (say, if the correlation coefficient does not meet a validation criterion), as shown in figure 3(b). This case will not be considered in this work, vectors missing due to vector post-processing will be filled with an interpolation algorithm. However, the approach that will be described in this section is agnostic to the nature of the missing vectors and should still handle these cases. Note figure 3 encodes the missing vectors as NaN (Not a Number) values, shown as a deep blue color in the contours. NaN values will encode missing vectors throughout this manuscript.

When vectors are missing due to masking (figure 3(a)), the integration of a pressure ray needs to be performed starting at the edges of the masked out region. Furthermore, the masked region can have concave regions where some of the rays will start and end without ever reaching the borders of the domain. Therefore, a more general approach to compute the weights $w_i^{\,j}$ considering arbitrarily masked border shapes is necessary to deal with realistic PIV datasets, especially in 3D volumes where the masked borders can have a large number of combinations of missing data topologies that would be intractable to handle manually in code.

Here a method will be introduced for the computation of the weights for arbitrary domain shapes. Only the 'face-crossing' scheme will be described in the body of this manuscript, for both 2D and 3D grids. The 'cell-centered' scheme is described in the supplementary material only for 2D grids, as the 3D grid integrals were exceedingly cumbersome to compute in the 'cell-centered' scheme with arbitrary boundaries.

2.2.1. Arbitrary boundary weights for 'face-crossing' scheme in 2D.

In order to compute the weights corresponding to interior cells with potentially arbitrary boundaries, let us consider the source field $\vec{f}$ has already been computed and it contains NaN values at some cells. If a cell has a NaN neighbor, rays coming from the NaN direction must be treated as boundaries and therefore must be assigned to $w_j^n$ (i.e. the nth iteration). In the face crossing scheme, a ray coming from the east (E) cell will have a path that passes through the east cell as well as one and only one of its neighbors (i.e. it must pass through either EN, EE or ES in figure 4(a)).

Zoom In Zoom Out Reset image size

We can define a boolean value b_j for each cell, defining whether there is missing data in the corresponding cell for the source field $\vec{f}$:

$$\begin{align} b_j = \begin{cases} 0 \;\; \textrm{if } \vec{f}\left(j\right) \neq \textrm{NaN} \\ 1 \;\; \textrm{if } \vec{f}\left(j\right) = \textrm{NaN}. \end{cases} \end{align} \tag{ 13 }$$

We also define $\tilde{b}_j = 1-b_j$ as the boolean NOT operation on b_j to define where data is available instead.

Thus, when considering all possible ray angles θ, we need to assess six possible cases depending on where the rays come from and whether the corresponding adjacent cells have missing data or not. The cases are enumerated in figure 4(b). For example, case A represents rays coming at an angle such that all rays come from the ES cell but some rays miss the C cell (and therefore must not be counted). Case B represents rays coming at an angle such that all rays coming from the ES cell also enter the C cell. Case C is the complementary case, where the rays at the same angle come from the EE cell and cross the C cell. Cases D, E and F are complementary cases to C, B and A, respectively, but involving the EN cell instead.

For each of the cases, a different mathematical expression representing the projection of the edge highlighted in blue can be defined. The total ray projection for the $(n+1)$th iteration then becomes:

$$\begin{align} \begin{aligned} A_E^{n+1} & = \tilde{b}_E \left[ \underbrace{\tilde{b}_{ES} \int_{-\pi/2}^\alpha \Delta y \text{cos}\theta \textrm{d}\theta}_{\textrm{Case A}} + \underbrace{\tilde{b}_{ES} \int_{-\alpha}^0 \Delta x \text{sin}|\theta| \textrm{d}\theta}_{\textrm{Case B}} + \underbrace{\tilde{b}_{EE} \int_{-\alpha}^0 \left(\Delta y \text{cos}\theta - \Delta x \text{sin} |\theta|\right) \textrm{d}\theta}_{\textrm{Case C}} \dots\right. \\ & \quad \left.+\ \underbrace{\tilde{b}_{EE} \int_{0}^\alpha \left(\Delta y \text{cos}\theta - \Delta x \text{sin} \theta\right) \textrm{d}\theta}_{\textrm{Case D}} +\underbrace{\tilde{b}_{EN} \int_{0}^\alpha \Delta x \text{sin}|\theta| \textrm{d}\theta}_{\textrm{Case E}} +\underbrace{\tilde{b}_{EN} \int_{\alpha}^{\pi/2} \Delta y \text{cos}\theta \textrm{d}\theta}_{\textrm{Case F}} \right] \end{aligned} \end{align} \tag{ 14 }$$

where $\alpha = \textrm{atan}(\Delta y/\Delta x)$. Note that $\tilde{b}_E$ multiplies the entire expression, to indicate that if there is missing data in the E cell then this entire weight evaluates to zero. When the data in the E cell is missing, then the current cell being evaluated (cell C) must be the starting point for all rays going westward and therefore is itself a boundary. Therefore, there must also be a boundary contribution from the cell C of the form:

$$\begin{align} A_C^n = 2\Delta y\left(b_E + b_W\right) + 2\Delta x\left(b_N+b_S\right). \end{align} \tag{ 15 }$$

This form considers all four cells surrounding cell C, and the constants $2\Delta x$ and $2\Delta y$ are calculated the same way as shown in equation (5). We can evaluate the integrals in equation (14) to obtain the expressions:

$$\begin{align} A_E^{n+1} = \tilde{b}_E \left[c_{xy} \left(\tilde{b}_{ES} + \tilde{b}_{EN}\right) + c_{xx} \tilde{b}_{EE}\right] \end{align} \tag{ 16 }$$

where the constants c_jk are defined as:

$$\begin{align} \begin{aligned}c_{xx} = 2\left(\Delta d - \Delta x\right) \\ c_{yy} = 2\left(\Delta d - \Delta y\right) \\ c_{xy} = \Delta x + \Delta y - \Delta d \end{aligned} \end{align} \tag{ 17 }$$

and $\Delta d = \sqrt{\Delta x^2 + \Delta y^2}$. These constants come from evaluating the integrals in equation (14).

Similarly, the ray integral corresponding to the boundary values for iteration n has the following form:

$$\begin{align} A_E^{n} = \tilde{b}_E \left[c_{xy} \left(b_{ES} + b_{EN}\right) + c_{xx} b_{EE}\right]. \end{align} \tag{ 18 }$$

To compute the weights corresponding to the west (W) cell, one simply swaps E and W in equations (16) and (18). For the north (N) and south (S) cells, the constant c_xx also needs to be swapped to c_yy :

$$\begin{align} A_N^{n+1} & = \tilde{b}_N \left[c_{xy} \left(\tilde{b}_{EN} + \tilde{b}_{WN}\right) + c_{yy} \tilde{b}_{NN}\right] \end{align} \tag{ 19 }$$

$$\begin{align} A_N^{n} & = \tilde{b}_N \left[c_{xy} \left(b_{EN} + b_{WN}\right) + c_{yy} b_{NN}\right]. \end{align} \tag{ 20 }$$

For a uniform rectangular grid, the constants c_xx , c_yy and c_xy can be trivially precomputed using equation (17). Then, the ray counting process to determine $A_j^i$ reduces to a conditional sum for each cell in the domain using the values b_j that define whether data for the function $\vec{f}$ is available in the corresponding neighboring cells. The weights $w_j^i$ then are found using equation (10).

2.2.2. Arbitrary Boundary Weights for 'Face-Crossing' Scheme in 3D.

In the three-dimensional case, the computation of weights is significantly more involved, and therefore full detail of the integrals computated will be provided in the supplementary material. In this manuscript, we describe the key steps involved in the weight computation in 3D. As the rays now come from all spherical angles, we have to compute double integrals with curved bounds. The complexity of the integrals is also increased because the projections of the neighboring faces onto the faces of the cell C are no longer lines, but triangles, trapezoids or rectangles depending on the values of the spherical angles.

The problem is still tractable, as will be demonstrated herein. Not all integrals will have a closed form, thus requiring them to be numerically computed. However, this computation only occurs once and only a handful of coefficients needs to be stored in memory. Similarly to the 2D case, this method has the advantage of treating all possible combinations of boundary conditions in 3D, which is combinatorially intractable. In the 2D 'face-crossing' scheme, there are $2^{12} = 4096$ possible combinations of boundary conditions. In 3D, the number of combinations is $2^{32}\approx 4$ billion. Evidently, most combinations in 3D would be degenerate cases, but it is still better to have a system that can treat any possible combination of boundary conditions than having to treat each boundary condition possibility within a pre-programmed case structure.

First, we will find the value of A_tot , considering all rays coming from all directions. For this, we will define a spherical coordinate system with polar angle θ and azimuthal angle φ, according to the right-handed coordinate system described in the ISO 80000-2:2019 convention. The definitions of the angles θ and φ are depicted in figure 5(a).

Zoom In Zoom Out Reset image size

Counting the total number of rays crossing the north face (N), which has a positive y normal, is straightforward. Defining a unit vector in the ray direction $\hat{r}$:

$$\begin{align} \hat{r} = \begin{bmatrix} \text{sin}{\theta} \text{cos}{\phi} \\ \text{sin}{\theta} \text{sin}{\phi} \\ \text{cos}{\theta} \end{bmatrix}. \end{align} \tag{ 21 }$$

The projection of the face N in the direction $\hat{r}$ is found by multiplying the surface area of N ($\Delta x \Delta z$) by the dot product $\hat{y}\cdot\hat{r}$. To find the total projection from all spherical angles A_N , one needs to integrate:

$$\begin{align} A_N = 4 \int_0^{\pi/2} \int_0^{\pi/2} \Delta x \Delta z \left(\hat{y} \cdot \hat{r}\right) \text{sin}{\theta} \textrm{d}\theta \textrm{d}\phi. \end{align} \tag{ 22 }$$

The factor of 4 is present due to the symmetry of the problem and the defined bounds of the angles θ and φ. Alternatively, one could make the upper bounds of both angles π and remove the factor of 4, yielding the same results.This integral can be evaluated analytically:

$$\begin{align} A_N = 4 \Delta x \Delta z \int_0^{\pi/2} \int_0^{\pi/2} \left(\text{sin}{\theta} \text{sin}{\phi} \right) \text{sin}{\theta} \textrm{d}\theta \textrm{d}\phi = \pi \Delta x \Delta z. \end{align} \tag{ 23 }$$

Given the symmetry of the problem, the total projection for the other directions ($A_E, A_W, A_F, A_B, A_S$) can be found by swapping the grid spacings Δx and Δz by the appropriate values. Doing this yields a total projection area of:

$$\begin{align} A_{tot} = 2\pi \left(\Delta x \Delta y +\Delta x \Delta z +\Delta y \Delta z \right). \end{align} \tag{ 24 }$$

Within the same theme of utilizing a boolean value b_j to define whether a cell j has NaN as a source term $\vec{f}$ (and therefore must be treated as a boundary), we can evaluate the area integrals for the interior and boundary contributions in 3D in a similar form as it was performed in 2D:

$$\begin{align} \begin{aligned} A_N^{n+1} & = \tilde{b}_N \left[c_{xy} \left(\tilde{b}_{EN} + \tilde{b}_{WN}\right) + c_{yz} \left(\tilde{b}_{NF} + \tilde{b}_{NB}\right) + c_{yy} \tilde{b}_{NN}\right] \\ A_N^{n} & = \tilde{b}_N \left[c_{xy} \left(b_{EN} + b_{WN}\right) + c_{yz} \left(b_{NF} + b_{NB}\right) + c_{yy} b_{NN}\right] \end{aligned}. \end{align} \tag{ 25 }$$

The pattern is also the same, where (N and S), (W and E) and (F and B) are swappable in the expressions above, and the corresponding constants c_jk must correspond to the appropriate directions (i.e. y for N and S, x for E and W, and z for F and B). This means a total of 6 constants c_jk must be calculated.

For the far faces (j = k), the area projection is displayed in figure 5(b), where P₁, P₂ are points on one of the edges of the far face and $P_1^{\prime}$, $P_2^{\prime}$ are the projections of these points on the N face through the $\hat{r}$ direction. Generalizing the projection for arbitrary coordinates $(j, k, l)$ results in the following integral that finds the contribution of the far face to the area projection A:

$$\begin{align} \begin{aligned} c_{jk} \left(j = k\right) & = \int_{\textrm{atan}\left(\Delta j/\Delta \ell_2\right)}^{\pi/2} \pi\left(\Delta \ell_2 \Delta \ell_{1,s} - \Delta j \Delta \ell_{1,c}\right) \\ & \quad + \dots 2\left[\Delta j\Delta \ell_{1,c} \;\textrm{atan}\frac{\Delta j}{\Delta \ell_{1,s}} - \Delta \ell_2\Delta \ell_{1,s} \;\textrm{atan}\frac{\Delta j}{\Delta \ell_{1,s}} \right] \textrm{d}\phi \end{aligned} \end{align} \tag{ 26 }$$

where Δj is the grid spacing in the j direction and $\Delta \ell_1$ and $\Delta \ell_2$ are the grid spacings in the remaining two directions. Furthermore, $\Delta \ell_{1,c} = \Delta \ell_1 \text{cos}{\phi}$ and $\Delta \ell_{1,s} = \Delta \ell_1 \text{sin}{\phi}$. Although $\Delta \ell_1$ and $\Delta \ell_2$ do not seem to be interchangeable in equation (26), the result from the integral is identical when $\Delta \ell_1$ and $\Delta \ell_2$ are swapped. Further details on how to achieve the integral form of equation (26) can be found in the supplementary material.

For the faces forming dihedrals (i.e. for c_jk where j ≠ k), the projections on the N face can form triangles or trapezoids, depending on the values of the angles θ and φ. There are three cases that must be considered, enumerated in figure 5(c). These cases result in four integrals:

$$\begin{align} & c_{jk} \left(j\neq k\right) = I_{A} + I_{B_1}+ I_{B_2}+ I_{C} \end{align} \tag{ 27 }$$

$$\begin{align} I_{A} & = \frac{\Delta \ell^2}{2} \textrm{atanh} \left( \frac{\Delta j^2 \Delta k^2}{2\Delta \ell^4 + 2\Delta \ell^2 \Delta j^2 +2\Delta \ell^2 \Delta k^2 + \Delta j^2 \Delta k^2}\right) \end{align} \tag{ 28 }$$

$$\begin{align} I_{B_1} & = \int_0^{\textrm{atan}\left(\Delta j/\Delta \ell\right)} \left\{\Delta \ell \Delta k_s\left( \textrm{atan} \frac{\Delta \ell}{\Delta k_c} - \frac{\Delta \ell \Delta k_c}{\Delta \ell^2 + \Delta k_c^2} \right) \right.\nonumber\\ & \left. \quad + \dots \Delta k_c \Delta k_s\left[\text{ln} \left(\frac{\Delta k_c}{\sqrt{\Delta k_c^2 + \Delta \ell^2}}\right) + \frac{\Delta \ell^2}{2\left(\Delta \ell^2 + \Delta k_c^2\right) }\right]\right\} \textrm{d}\phi \end{align} \tag{ 29 }$$

$$\begin{align} I_{B_2} & = \int_{\textrm{atan}\left(\Delta j/\Delta \ell\right)}^{\pi/2} \left\{\Delta \ell \Delta k_s\left( \textrm{atan} \frac{\Delta j}{\Delta k_s} - \frac{\Delta j \Delta k_s}{\Delta j^2 + \Delta k_s^2} \right)\right. \nonumber\\ & \left. \quad + \dots \Delta k_c \Delta k_s\left[\text{ln} \left(\frac{\Delta k_s}{\sqrt{\Delta k_s^2 + \Delta j^2}}\right) + \frac{\Delta j^2}{2\left(\Delta j^2 + \Delta k_s^2\right)}\right]\right\} \textrm{d}\phi \end{align} \tag{ 30 }$$

$$\begin{align} I_{C} & = \frac{1}{2} \int_{\textrm{atan}\left(\Delta j/\Delta \ell\right)}^{\pi/2} \frac{\Delta j \Delta k^2 \text{sin}{\phi} \left(2\Delta \ell \text{sin}{\phi} - \Delta j \text{cos}{\phi}\right)}{\Delta j^2 + \Delta k^2 \text{sin}^2 \phi} \textrm{d}\phi \end{align} \tag{ 31 }$$

where $\Delta \ell$ represents the grid spacing in the direction not considered ($\ell \neq \{j,k\}$) and $\Delta k_c = \Delta k \text{cos}{\phi}$ and $\Delta k_s = \Delta k \text{sin}{\phi}$.

Note equation (28) has a closed form, whereas the remaining integrals need to be computed numerically. The forms provided are the simplest analytical expressions the authors were able to derive. Although not straightforwardly evident at a first glance, $c_{jk} = c_{kj}$. This is because the projection of face j onto face k for any viewing angle always has the same area as the projection of face k onto face j. Once again, further details on the computation of the integrals can be found in the supplementary material.

Although the integrals in equations (26) and (29)–(31) and must be numerically computed, this computation only needs to be performed once prior to running the matrix omnidirectional solver. Computing all integrals to $1\times10^{-7}$ precision takes about 3 s for arbitrary Δx, Δy and Δz, which was deemed a reasonable investment for solving an arbitrarily-spaced fully three-dimensional grid. For common ratios of Δx, Δy and Δz, these integrals can also be precomputed and stored. A special case is shown in the next section.

2.2.3. 3D 'Face-Crossing' Constants for Isotropic Grids.

For the special case where the grid spacing is equal in all directions (i.e. $\Delta x = \Delta y = \Delta z = \Delta$), we can perform further simplification of the integrals and attain a closed form for c_jk , which then can be added to the solver code as a fixed constant for isotropic grids. Below we can find the values of these constants:

$$\begin{align} c_{jk} \left(j = k\right) & = \Delta^2 \left\{\text{ln} \left(\frac{4}{3} \right) - \sqrt{2} \; \textrm{atan} \left(\sqrt{2}\right) \right. \nonumber\\ & \left. \quad +\ \frac{\sqrt{2}}{2} \textrm{atan} \left(\frac{10\sqrt{2}}{23} \right) + \pi \left(\sqrt{2} -1\right) \right\} \nonumber\\& \qquad\qquad\qquad\qquad\qquad \textrm{(Isotropic Grid)} \end{align} \tag{ 32 }$$

$$\begin{align} c_{jk} \left(j\neq k\right) & = \Delta ^2 \left\{\frac{\pi}{2} - \frac{1}{4}\text{ln} \left(\frac{4}{3}\right) \right. \nonumber\\ & \left. \quad -\ \frac{\sqrt{2}}{4} \textrm{atan}\left(\frac{10 \sqrt{2}}{23}\right)- \frac{\sqrt{2}}{2} \; \textrm{atan}\left(\sqrt{2}\right)\right\}.\nonumber\\ & \qquad\qquad\qquad\qquad\qquad \textrm{(Isotropic Grid)} \end{align} \tag{ 33 }$$

Numerically, these constants have the following values:

$$\begin{align} c_{jk} \left(j = k\right) & = 0.6277684243304 \Delta ^2\; \; \textrm{(Isotropic Grid)} \end{align} \tag{ 34 }$$

$$\begin{align} c_{jk} \left(j\neq k\right) & = 0.6284560573148 \Delta ^2\; \; \textrm{(Isotropic Grid)} \end{align} \tag{ 35 }$$

First, a decaying Taylor vortex is considered, which was previously used by Charonko et al (2010) to assess error propagation in pressure integration algorithms. The Taylor vortex has an analytical solution for pressure that enables the definition of a ground truth. The velocity field of the decaying Taylor vortex centered at the origin is defined in cylindrical coordinates as:

$$\begin{align} u_\theta \left(r\right) = \frac{H}{8\pi} \frac{r}{\nu t^2} \text{exp} \left(-\frac{r^2}{4\nu t}\right) \end{align} \tag{ 36 }$$

where u_θ is the tangential velocity, r is the radial position and H is a measure of the angular momentum in the vortex. The pressure field for this vortex can be found by solving the Navier–Stokes equations:

$$\begin{align} P \left(r\right) = -\frac{\rho H^2}{64\pi^2 \nu t^3} \text{exp} \left(-\frac{r^2}{2\nu t}\right) + P_\infty \end{align} \tag{ 37 }$$

where $P_\infty$ is an integration constant and is defined as the far-field pressure.

For this study, the parameters $H = 10^{-6}$ $\textrm{m}^2$, $\nu = 10^{-6}$ $\textrm{m}^2/\textrm{s}$, $\rho = 1,000$ $\textrm{kg}/\textrm{m}^3$ were used. The characteristic scales for length ($L_0 = \sqrt{H}$), velocity ($U_0 = \nu/\sqrt{H}$), time ($T_0 = H/\nu$) and pressure ($P_0 = \rho U_0^2$) are also defined to normalize the values prior to sending the data to the solvers.

In order to generate velocity fields $\vec{u}_{noisy}$ contaminated with spatially uncorrelated noise, a noisy vector field $\vec{n}$ with unit magnitude and random angle is sampled from an uniform angular distribution and scaled by a factor s. It is then multiplied pointwise by the magnitude of the instantaneous velocity vectors $||\vec{u}||$, following a previous work (Charonko et al 2010):

$$\begin{align} \vec{u}_{noisy} = \vec{u} + s\vec{n} ||\vec{u}||. \end{align} \tag{ 38 }$$

The scaling factor s is then varied between 0 and 0.1 (0% and 10% noise intensity) for this analysis. In all cases examined, a POD-based denoising scheme retaining only the most energetic modes in the time series was used to reduce noise in the velocity fields, following the previous work by Charonko et al (2010). Modes with less than 1% of the total energy were discarded, resulting in approximately 6 modes retained for pressure reconstruction.

Five pressure reconstruction algorithms were tested in this part of the work, namely:

• (i)  
  Pressure-Poisson solver with Neumann boundaries;
• (ii)  
  Omnidirectional integration, 'virtual boundary' (Liu and Katz 2006);
• (iii)  
  Omnidirectional integration, 'rotating parallel ray' (Liu et al 2016);
• (iv)  
  Matrix-based omnidirectional solver, 'face-crossing' scheme;
• (v)  
  Matrix-based omnidirectional solver, 'cell-centered' scheme;

For the Pressure-Poisson case (i), we defined an additional Dirichlet boundary condition at one single mesh point at a corner of the domain to ensure the problem is well-defined. To define the source terms, the 'standard' approach as defined in Charonko et al (2010) was used. This method used the full Navier–Stokes equations for Neumann boundaries, and direct calculation of the divergence of each term for the source term, preserving the time-derivatives and viscous terms with no application of an incompressible flow assumption. This was chosen for consistency with the forms used in the omni-directional solvers, but as shown in Charonko et al (2010) and Nie et al (2022), alternate forms that reduce numerical error propagation can sometimes result in better performance.

For the 'virtual boundary' case (ii), a total of $2 N_x N_y$ rays were used in the boundary, whereas for the 'rotating parallel ray' approach (iii) a ray spacing $\Delta s/\Delta x = 0.4$ was used and Δθ was four times the value provided in formula (73) in Liu and Moreto (2020).

To compare computational efficiency and accuracy of the algorithms tested, a 2D version of the algorithms described is implemented on a square domain of variable number of grid points containing the Taylor vortex flow field. The omnidirectional integration algorithms (ii) and (iii) are implemented in Matlab using the parallel computing toolbox, and 8 nodes of an Intel Xeon Gold 5218R CPU (2.10 GHz) were used to accelerate the computations. The algorithms requiring a matrix solver (i), (iv), (v) employed Matlab's default linear equation system solver (mldivide or 'backslash' operator). Finally, in the error tests considered in this section, a total of 26 Taylor vortex snapshots at times spanning $0.05\lt t/T_0\lt0.3$ using a 101² grid that spans from $-3\lt x/L_0\lt3$ and $-3\lt y/L_0\lt3$ were employed.

A summary of the RMS pressure errors, as a percentage of the maximum absolute pressure in the field, are presented in figure 6(a) for all solvers considered. It becomes quickly evident that the pressure-Poisson solver presents very poor accuracy, providing incorrect results with an RMS error of ${\sim}5\%$ of the peak suction even when the underlying velocity field contains no noise. The error of the pressures obtained by the pressure-Poisson solver then shoots up very quickly as a function of velocity error. The two omnidirectional methods ('virtual boundary' and 'parallel ray') perform better than the pressure Poisson method, providing a somewhat linear relationship with the velocity error with a zero crossing at approximately zero pressure error. However, we can see the statistics from both methods fluctuate for different velocity errors, due to the limited number of rays employed during the integration. We did not run more cases for converged error statistics for the traditional omnidirectional methods, as they are highly time intensive. The matrix omnidirectional methods proposed in this manuscript, on the other hand, offer the lowest errors; the 'cell-centered' approach performing slightly better. Furthermore, the statistics from the matrix methods are better converged, as the continuous integrals utilized mimic the limit of infinite rays.

Zoom In Zoom Out Reset image size

When the computational cost to evaluate pressure from the velocity fields is considered, there are several points that will be further elaborated throughout this manuscript. For this first analysis on the 2D Taylor vortex, we consider the computational time required to compute one iteration of each one of the five methods examined in figure 6(b) for a variety of different grid sizes ranging from 11² to 2001² in 2D. We also include in figure 6(b) the computational times for the GPU implementation in 2D and 3D (for grids of 51³ to 501³), shown as dot symbols. The only method that requires a single iteration to converge is the pressure-Poisson method, and all other methods will require multiple iterations to converge. Less than 30 iterations were required for the Taylor vortex problem examined in this section for all omnidirectional methods, but further discussion is presented in section 3.4 regarding the convergence behavior for more complex cases, where hundreds of iterations may be required to achieve satisfactory convergence. In our numerical experiments, however, the matrix omnidirectional methods always converge in less iterations to a given residual than the traditional omnidirectional methods.

Figure 6(b) demonstrates how significant the computational time reduction is for the matrix omnidirectional methods discussed herein, requiring a logarithmic scale to fit the times from the different methods. One iteration of the matrix methods requires slightly more time to solve in Matlab than the pressure-Poisson equation, but the times are of the same order. On the other hand, the traditional omnidirectional integration schemes require about 20 times longer per iteration for a given grid size in 2D for the discretization chosen here. The number of iterations required to convergence of the pressure to a certain criterion is at least the same as the matrix methods proposed herein.

When we consider our GPU implementation of the 2D matrix omnidirectional solver, we observe a ∼3–3.5× speedup in comparison to our Matlab implementation. For 3D grids, we only implemented the matrix omnidirectional solver in the GPU, and curiously we find that the 3D grids solve significantly faster for a given precision in the CG solver for the same number of grid points, averaging about ∼9× faster than a 2D grid with the same total number of grid points. We have not explored the underlying reason for this drop in computational complexity, but it is likely related to the number of iterations required for the CG solver to propagate information through the domain, which is of the order of the largest dimension in the domain. However, this effect is particularly exciting because it means that a speedup of approximately 3 orders of magnitude (∼500–1000× faster) can be achieved when compared to 2D grids with the traditional omnidirectional methods, which will enable the tackling of large 3D grids that are becoming highly relevant to the experimental fluids community.

Although a 3D implementation of the traditional omnidirectional integration was not currently examined in this work, an important time complexity scaling is evident by examination of the work presented by Wang et al (2019): Using $\mathcal{N}$ to denominate the largest size of the grid between the x, y and z directions (assuming a box-shaped domain), the number of rays cast is proportional to the projection area of the domain box analyzed at a given spherical angle (i.e. ${\sim}\mathcal{N}^2$). The number of operations required to compute each integral is ${\sim}\mathcal{N}$. Furthermore, the number of spherical angles required to tesselate the enclosing sphere to a given grid spacing is proportional to the area of the enclosing sphere (${\sim}\mathcal{N}^2$). Therefore, the number of computations required per iteration must be $t\sim \mathcal{O}(\mathcal{N}^5)$, which is a rather strong scaling. We find in our empirical tests of our 3D GPU implementation of the matrix omnidirectional method presented herein that the time complexity scaling is $t\sim \mathcal{O}(\mathcal{N}^{4})$, as shown in figure 6(b), which is significantly more favorable for larger grids. Note that the total number of points in the 3D grid scales as ${\sim}\mathcal{N}^{3}$.

A more direct comparison with the performance described by Wang et al (2019) can be performed by comparing the performance of the graphics cards used in both studies. The RTX A4000 used in our work is capable of 19.2 TFLOPS at single-precision, whereas the RTX 2080 Ti used by Wang et al (2019) is capable of 13.45 TFLOPS. We performed our calculations in double-precision, which is reported by Wang et al (2019) to take approximately 3 min per iteration for a $100\times 47\times 38$ grid in their machine. A cutout of the JHU turbulence data described in section 3.2 was generated with the same $100\times 47\times 38$ grid size, and one iteration of the omnidirectional matrix 'face-crossing' scheme takes on average 0.082 s. Therefore, the GPU implementation of our method for 3D grids is estimated to be about 1500× faster than the state of the art (Wang et al 2019) at this specific grid size and likely even faster for larger grid sizes due to the more favorable time complexity scaling of $t\sim \mathcal{O}(\mathcal{N}^{4})$.

In this section, the performance of the omnidirectional matrix method will be assessed for a more complex flow field for which the exact values of pressure and velocity are known. The 'forced isotropic turbulence' direct numerical simulation (DNS) dataset provided by the Johns Hopkins University (Li et al 2008) was chosen for this analysis. More specifically, the 'isotropic1024coarse' database was sampled at the 40th frame (time = 0.2535) and adjacent frames for the computation of the time derivatives.

Grid resolutions of 50³, 100³, 200³ and 500³ were considered in this study. The data cutouts taken are not decimated or downsampled to the specified grid resolutions considered to minimize truncation error in the derivative computations. Instead, the cutouts have the same grid spacing as the original data set and therefore the 50³ box is truncated from a corner of the original 1024³ set and is four times smaller in physical size than the 200³ box. This also enables the demonstration of how the omnidirectional matrix method handles incomplete data boundaries, since none of the cases analyzed comprise of periodic boundaries. As boundary conditions are not defined in the omnidirectional matrix method, this poses no additional implementation challenge.

While sampling the JHU turbulence dataset, we noticed that the momentum equations are not balanced when computing the finite difference gradients of the available velocity/pressure fields in an Eulerian grid, even when including the forcing terms; for accuracy levels from second-order to sixth-order. This slight imbalance appears to be related to the fact the original database was solved with a spectral method. We believe that this effect acts as an unnecessary confounding factor that can contaminate the obtained pressure fields but does not pertain to the accuracy of the pressure solver we propose, but to missing information in the underlying velocity field. Therefore, in order to generate noise-contaminated velocity fields based on the correct underlying pressure field, we perform the following operations:

• (i)  
  Generate the uncontaminated source field $\vec{f}_0$ based on the true pressure with central finite differences: $\vec{f}_0 = \nabla P_{truth}$;
• (ii)  
  Compute the source term based on the convective derivative $\vec{f}_u = -\rho (D\vec{u}/Dt)$ based on the velocity fields using second-order accurate central finite differences;
• (iii)  
  Compute a residual source term $\vec{r} = \vec{f}_0-\vec{f}_u$;
• (iv)  
  Add noise to the velocity field ($\vec{u}_{noisy}$) and compute a noise-contaminated source term $\vec{f}_{u,n} = -\rho (D\vec{u}_{noisy}/Dt)$;
• (v)  
  Correct the source term with the residual: $\vec{f}_{0,n} = \vec{f}_{u,n} + \vec{r}$;
• (vi)  
  Integrate the noise-contaminated source field $\vec{f}_{0,n}$ to obtain pressure.

With this process, the source field $\vec{f}_{0,n}$ computed from the velocity field contaminated with noise of zero magnitude matches $\nabla P_{truth}$ and has no error related to the approximated derivatives. The noise-contaminated velocity fields $\vec{u}_{noisy}$ are generated differently than the previous section, as the fields are three dimensional. The equation below adds multiplicative noise to the velocity fields scaled by a factor s, between 0 and 0.1 (0% and 10% noise intensity). Note the noise added in equation (39) is slightly different than equation (38) in that it only rescales the vectors without adding a component perpendicular to them. Due to the usage of a single field to estimate pressure, no noise reduction scheme (such as the POD-denoising used in the previous section) was employed.

$$\begin{align} \vec{u}_{noisy} = \vec{u} \left(1+s\vec{n}\right). \end{align} \tag{ 39 }$$

Figure 7 displays the contours of pressure for two turbulence box sizes, 200³ and 500³, which were iterated to a 10⁻³ momentum residual as defined in equation (41) and further discussed in section 3.4. In each set of cases, the slice is taken at the middle of the respective domain. We can compare the true pressure in figures 7(a) and (f) against the results obtained with the matrix omnidirectional method at 0% noise (b), (g) and 5% noise (d), (i). The error fields (truth minus 'matrix-omni' result) are shown in figures 7(c), (e), (h) and (j) in a tighter color scale. It is evident from the contours in figure 7 that the matrix omnidirectional results closely approximate the true pressure fields in both cases even with the presence of fairly large velocity noise. A velocity field noise with a 5% standard deviation is representative of the quality that might be achieved with particle image velocimetry for the more challenging experiments where illumination intensity is very low. Nevertheless, we still observe the shape and magnitude of the contours in the noisy cases is close to the true pressure, up to quite significant detail.

Zoom In Zoom Out Reset image size

Several cases with varying velocity noise were tested with this data set and the RMS pressure error is plotted in figure 8 as a function of velocity noise for various grid sizes. We note a similar linear trend as a function of velocity noise as observed in the 2D Taylor vortex shown in section 3.1, however the RMS pressure error is smaller in the JHU turbulence case. This may be related to the more uniformly distributed magnitudes of the velocity vectors in the JHU database, which reduces the error of the velocity gradients required to compute the source terms $\vec{f}$. The generally low error achieved for different grid sizes demonstrates the ability of the matrix omnidirectional method proposed herein to provide accurate pressure estimates from PIV velocity fields even for very large 3D mesh sizes.

Zoom In Zoom Out Reset image size

In previous sections, we argued that the matrix omnidirectional method we propose simply takes the 'rotating parallel ray' approach to the limit of infinite rays. From the integral formulation proposed in section 2.1, we define the matrix equation weights considering the infinite number of rays coming from all possible angles θ between 0 and 2π and all possible ray shifts crossing any given cell in the domain. Therefore, the results from the computations using the 'rotating parallel ray' approach with increasingly smaller spacing between rays and increasingly more angles should approach the result from the matrix method proposed herein.

This behavior is evidenced in figure 9, where the 'rotating parallel ray' algorithm is run on the 2D Taylor vortex data with 1% velocity error on a 101² grid for 30 iterations of the Matlab implementation. The blue solid line, showing the RMS pressure error as a percentage of the peak suction at the core of the vortex, evidences this convergence behavior as the number of rays (integration paths) is increased. However, note the logarithmic x-scale, demonstrating a large number of integration paths is required to really achieve this convergence. Evidently, this comes at a high computational cost, which is evidenced by the red solid line in figure 9. Compared to the time required to run the 30 iterations of the matrix omnidirectional method (shown as a red dashed line, ${\sim}2.07$ s), the cost to run the 'rotating parallel ray' algorithm to the same accuracy as the matrix omnidirectional method is very high, of order 10 min, even for a small 101² domain.

Zoom In Zoom Out Reset image size

As discussed in section 2, equation (12) is an iterative update equation for the pressure field. It updates both boundary nodes and internal nodes for the $(n+1)$th iteration. Here we discuss the observed convergence properties of this update method.

In order to define convergence, two metrics were examined. The pressure update metric is defined as:

$$\begin{align} \varepsilon_P^{n+1} = \frac{||\left\{P^{n+1}\right\} - \left\{P^{n}\right\}||_2}{||\left\{P^{n+1}\right\}||_2} \end{align} \tag{ 40 }$$

where $||\cdot||_2$ is the L2 norm of the field considered. Alternatively, the residual of equation (12) (momentum residual) can be used as a convergence metric:

$$\begin{align} \varepsilon_R^{n+1} = \frac{||\left[W^{\,n+1}\right] \left\{P^{n+1}\right\} - \left[W^{\,n}\right] \left\{P^{n+1}\right\} -\left\{S\right\}||_2}{||\left\{S\right\}||_2} \end{align} \tag{ 41 }$$

where the denominator using only the source term assumes the first residual of equation (12) was computed with a zero initial guess for $\{P^{n = 0}\} = 0$. The momentum residual of equation (41) is readily available in the CG solver implementation, as it needs to be calculated when the CG solver is initialized.

Although both residuals could be used for assessing convergence, the momentum residual ε_R was found to be a slightly more conservative convergence criterion (i.e. its decay rate is the same as ε_P but its magnitude is slightly larger) and therefore it will be the quantity used for all discussions through this manuscript.

Due to the high computational cost of the traditional omnidirectional integration method, previous references typically use a small number of iterations of the algorithm. For example, McClure and Yarusevych (2017) used 5 iterations in their examination of a 2D slice of a DNS cylinder flow database, whereas Liu and Katz (2006) claim to have observed convergence in less than 10 iterations in a 2D PIV measurement of a cavity shear flow. In the 3D implementation of the omnidirectional integration method by Wang et al (2019), only 3-4 iterations were used, which was deemed to have reached convergence by the authors. With the implementation proposed herein, the effect of the number of iterations on convergence can be examined for a very large number of iterations, as the computational cost per iteration is small even for large grids.

3.4.1. Effect of pressure magnitude at the boundaries.

First, we examine the convergence behavior for the 2D decaying Taylor vortex at a specified time ($t/T_0 = 0.1$) and varying the magnitude of the true pressure values at the boundary by changing the extents of the domain considered under the same grid resolution of 101² and no noise contamination of the underlying source fields.

The 'pressure magnitude' here is defined as the norm of the true pressures at the boundary $||P||_b$. It stands to reason that the further the final pressures at the boundary are from the initial guess, the longer it will take for the algorithm to update the boundary values to within a given residual. As the initialization of the pressure term $\{P^{n = 0}\} = 0$ has zero norm, this effect can be captured by the norm of the boundary pressures $||P||_b$, which is evidenced in figure 10. If the domain sampled from the Taylor vortex is very large, then the pressures at the boundary are close to zero (i.e. their norm $||P||_b\rightarrow 0$). Conversely, with a smaller windowing of the domain, the pressures at the boundary will have values further from zero, and their norm grows larger. As evidenced in figure 10(a), increasing the norm of the target boundary pressure slows down the convergence of the algorithm. The trend is clearer when the convergence criterion is fixed, as shown in figure 10(b), where ε_R was fixed to 10⁻⁴ and the number of iterations to convergence was counted for several values of $||P||_b$. The number of iterations to convergence appears to be a logarithmic function of the boundary norm $||P||_b$.

Zoom In Zoom Out Reset image size

This effect may explain why some of the former references observed rapid convergence in the traditional omnidirectional integration method. If the pressure field boundaries include a large region with pressures close to the ambient (or 'far field') pressure, then a fast convergence rate is attained because the boundary values do not need to be significantly updated. This also serves as a guideline for experimentation; i.e. the convergence rate will be faster when the far field is included in the velocity field measurements. However, this is not possible when the domain includes walls where a pressure distribution is present, as will be demonstrated with real data in section 4.

3.4.2. Grid size.

Now considering fields with complex pressure boundaries, we evaluate the effect of grid size on the number of iterations required to attain convergence. In this section, the same JHU isotropic turbulence data set examined in section 3.2 is considered with different grid sizes. The convergence behavior is summarized in figure 11. In the 2D cases shown in figure 11(a), the 2D slices are taken from the center of the 3D box of same size. It is evidenced in both figures 11(a) and (b) that the convergence behavior is somewhat complex with an initial polynomial phase and a rapid exponential convergence at the latter iterations. The 'knee' of the curves where the exponential phase begins seems to be different for each case analyzed, resulting in crossing of the convergence lines from each case. Generally, however, the convergence is faster for the smaller grids at the smaller residuals ε_R and less than 10³ iterations is required for most cases considered to converge to $\varepsilon_R\lt10^{-4}$.

Zoom In Zoom Out Reset image size

Curiously, the number of iterations required to compute the 3D cases (which include far more boundary points) is comparable to the 2D cases. The wall clock time required for the 2D cases, however, is far smaller; taking 310 s to run 1000 iterations for the 400² grid and 2863 s to run 1000 iterations for the 200³ grid using the GPU implementation of this algorithm.

3.4.3. Relationship Between Residual and Error.

When the source field $\{S\}$ used as an input to the omnidirectional matrix algorithm is contaminated by noise, it may be unnecessary to fully converge the pressure integration. Therefore, it is worthwhile to examine the relationship between the residual convergence and the error in the pressure fields for noisy cases to better understand where satisfactory convergence is attained.

Figure 12(a) presents the residual convergence for the 100³ grid case from the JHU dataset for three different noise magnitude cases: 0% noise, 2% noise and 5% noise. The first observation is that the cases with larger noise converge to lower residuals ε_R earlier than the uncontaminated case. Figure 12(b) demonstrates that the addition of noise in the underlying velocity fields makes the more converged cases plateau to a given RMS error in pressure when compared to the true pressure. For example, although $\varepsilon_R = 10^{-5}$ results in an RMS pressure error of 0.15% for the uncontaminated case, a $\varepsilon_R\sim10^{-4}$ is sufficient to attain the best pressure error of ${\sim}1\%$ for the 2% noise case and $\varepsilon_R\sim10^{-3}$ is sufficient to attain the RMS error of ${\sim}2\%$ for the 5% noise case. Note that under these conditions, somewhere between 100 and 500 iterations would already be sufficient for convergence, as shown in figure 12(a), and therefore a momentum residual of $\varepsilon_R\sim10^{-3}$ to 10⁻⁴ may be an appropriate stopping criterion for this algorithm under practical conditions where noise is present in the underlying velocity fields.

Zoom In Zoom Out Reset image size

3.4.4. Wall-clock time convergence.

It is also worth mentioning that the use of the Conjugate Gradient algorithm to solve the matrix equation, as implemented in the GPU version of this algorithm, presents the added advantage of enabling the use of the previous iteration's pressure field as a first guess for the CG algorithm. This significantly speeds up the solution time at the latter iterations of the algorithm, as is demonstrated in figure 13 for the 100³ grid from the JHU isotropic turbulence dataset. As it can be seen, after iteration ${\sim}300$ the time per iteration slowly decreases, reaching a negligible baseline of ${\sim}6$ milliseconds after iteration 1500 (corresponding to the various kernel launch overheads in the GPU). The peaks in the iteration times (blue lines) are due to periodic saving of the intermediary fields every 50 iterations.

Zoom In Zoom Out Reset image size

First, we consider the bluff body wake of a cylinder with a slanted afterbody examined by Zigunov et al (2020) at the slant angle of $\phi = 45^\circ$ and diameter-based Reynolds number of $Re_D = 25,000$. At these conditions, this flow presents an interesting hysteresis behavior, where the same Reynolds number supports two distinct wake flow topologies. In the 'vortex state', the flow separates at the leading edge of the slanted surface and reattaches on the slanted surface further downstream, forming a separation bubble and a counter-rotating vortex pair. In the 'wake state', the flow separates at the same point but does not reattach, forming a more traditional fully-separated bluff body wake. The 'wake state' is the default state for this wake at low Reynolds numbers, but if the Reynolds number exceeds a threshold value (of about 33 000), then the 'vortex state' forms and persists if the free stream velocity is then reduced (Zigunov et al 2022), allowing a range of Reynolds numbers to display two different flow topologies at the same Re_D .

The surface pressure distribution was measured using an Omega PX653 pressure transducer with a range of 0-0.1 inH2O, resulting in a C_p uncertainty of ±0.03 for the conditions analyzed. Five hundred planar PIV velocity fields were acquired at the center plane (streamwise-vertical) at a rate of 10 Hz, enabling the approximation of the Reynolds stresses required to compute the average pressure field.

The average pressure field is calculated by using the Reynolds-Averaged Navier Stokes equations, dropping the viscous terms:

$$\begin{align} f_i = \delta_{ij}\frac{\partial \bar{P}}{\partial x_j} = -\rho \left(\bar{u}_j\frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \left(\overline{u_i^{\prime} u_j^{\prime}}\right)}{\partial x_j}\right). \end{align} \tag{ 42 }$$

Since the underlying velocity fields were from planar PIV, the out-of-plane derivatives and velocity components were considered equal to zero. As the plane chosen is the symmetry plane of the flow, this is true for the mean quantities but not for the fluctuations.

The forcing vector field f_i is computed using equation (42) and integrated using the matrix omnidirectional method using the 'face-crossing' scheme with our GPU implementation. The PIV vector fields are 427×360 (∼154 k grid points) in size and masks were applied to remove regions without valid vectors (∼87 k grid points contain valid data). The first iterations took approximately 0.4 s, but about 3000 iterations were required to converge the boundaries; taking about 16 min to converge the pressure field with the matrix method.

Although the number of iterations and convergence time was large for this case, the results presented in figure 14 for the 'vortex regime' demonstrate the motivation for attaining convergence of the boundaries in the matrix omnidirectional method. The contours of pressure obtained with the matrix omnidirectional method in figure 14(a) are clearly a different field than the results from the pressure-Poisson method in figure 14(b), which was obtained by setting the boundary conditions to Neumann boundaries $(\nabla P \cdot \hat{n} = \vec{f}\cdot \hat{n})$. In both cases, the pressure at the top-right corner of the field, which is the closest to a free stream value, is set to $C_p = 0$.

Zoom In Zoom Out Reset image size

In the comparison with the transducer measurements shown in figure 14(c), it becomes clear that the matrix omnidirectional method presents a better match without any a priori knowledge (i.e. the only information provided was the $C_p = 0$ value at the top-right corner after pressure was integrated). The magnitude of the surface pressures is off by approximately $\Delta C_p = 0.3-0.5$ in the pressure-Poisson method. Furthermore, the low-pressure region towards the trailing edge $z/D\approx 1$ has incorrect topology in the pressure-Poisson solution of figure 15(b), incorrectly assigning a suction zone to the trailing edge flow. These aspects demonstrate the limitations of using the pressure-Poisson equation to obtain pressure estimates from PIV fields. Although the pressure-Poisson solution may seem adequate at a first glance without any other reference pressure field for comparison, it may lead to misleading interpretations about the physics of a problem. Therefore, the computational time investment for solving the omnidirectional matrix equation through a large number of iterations is likely justified if more accurate and robust pressure estimates are to be extracted from experimental PIV fields.

Zoom In Zoom Out Reset image size

For the sake of completeness, the pressure field for the 'wake regime' of this problem is also provided in figure 15(a), only for the matrix omnidirectional method. Note the remarkable agreement with the surface measurements in figure 15(b), likely due to the lower surface-normal pressure gradients. In other words, the pressure is not varying as much away from the surface as in the 'vortex regime', minimizing the effect of the missing vectors near the surface (due to laser reflections) and enabling a fairer comparison between the surface pressure from transducers and near-surface pressure from PIV. Note the pressure field topology in figure 15(a) matches the expected distribution for a fully separated bluff body wake, with the entirety of the wake presenting a moderate suction and the pressure recovering to free stream values towards the trailing edge.

In this application we demonstrate the possibility of applying this technique to obtain a 3D average pressure field for a compressible flow at $M_\infty = 2$ over a double-fin (data available in [Zigunov and Seckin, 2024]). The PIV fields are obtained using the newly-proposed 'scanning stereoscopic PIV' (SSPIV) technique (Zigunov et al 2023). This large volume ($53\times 66\times 100$ mm) with the $\lt$1 µm particles required for this supersonic flow is currently outside of the typical range for a tomographic PIV setup, which motivates the choice of this technique. In brief, the SSPIV technique comprises of a continuously scanning stereoscopic PIV system where the entire acquisition optics are slowly scanned throughout the interrogation volume. A stack of about 800 planes is then acquired over the course of a run, and each plane is assigned to a z location according to the timing of the laser Q-switch in relationship to a simultaneously-measured linear potentiometer attached to the traverse mechanism. This allows for the calculation of the average velocities and approximation of the Reynolds stresses by performing a spatio-temporal average, where PIV planes from z locations within a threshold (±1 mm) of a given z coordinate are averaged. Further details about the acquisition technique are provided in Zigunov et al (2023), and further details on the flow and analysis of the 3D flow structures observed can be found in Seckin et al (2023b). As a final note, we point out that the grid spacing is different in the in-plane directions and the out-of-plane direction, which makes the computation of the rectangular grid coefficients provided in section 2.2.2 particularly useful.

Before discussing the implementation details, we note that 30 planes are used for each spatiotemporal average plane, and to attain better statistical convergence of the averages and fluctuations the vectors on a $5\times 5\times 30$ grid centered at every grid point are used. More planes could not be acquired due to the limited blowdown time of the facility. This number of vectors is sufficient to converge the mean quantities, as discussed by Zigunov et al (2023), but the Reynolds stresses are likely not converged. The contribution of the Reynolds stresses in this flow field is small, and a comparison between the pressure fields with and without including the Reynolds stresses was performed and negligible difference is observed. For this demonstration, the results are reported with the Reynolds stress terms included in the computation of the source field.

To obtain the 3D average pressure field measurements from PIV we use the compressible RANS equations, as described by van Oudheusden (2008). This requires three key, but reasonable assumptions: (1) adiabatic flow, (2) ideal gas and (3) negligible effect of density fluctuations. The latter assumption (3) is reasonable for moderate Mach numbers. The ideal gas assumption (2) is also valid for this flow, as the minimum compressibility factor is estimated to be $PV/RT = 0.998$ (Su 1946). Given assumptions (1) and (2), we can find the temperature field T from the velocity field by knowing the total temperature $T_\infty$: ($T_\infty = 173$ K from isentropic relations)

$$\begin{align} \frac{T}{T_\infty} = 1+\frac{\gamma-1}{2} M_\infty^2\left(1-\frac{V^2}{V_\infty^2}\right) \end{align} \tag{ 43 }$$

where γ is the specific heat ratio and V is the velocity magnitude. Once the temperature field is estimated, the following conservative form of the compressible RANS momentum equations is used: (van Oudheusden 2008)

$$\begin{align} \left(\delta_{ij}+\frac{\overline{u_i u_j}}{RT}\right) \frac{\partial \text{ln} \left(\bar{P}/P_\infty\right)}{\partial x_j} = -\frac{1}{RT} \left( \frac{\partial \overline{u_i u_j}}{\partial x_j} - \frac{\overline{u_i u_j}}{T}\frac{\partial T}{\partial x_j}\right) \end{align} \tag{ 44 }$$

where R is the ideal gas constant and $\overline{u_i u_j} = \bar{u}_i\bar{u}_j + \overline{u_i^{\prime} u_j^{\prime}}$. Equation (44) can be integrated using the matrix omnidirectional method to obtain $\text{ln} (\bar{P}/P_\infty)$, and therefore $\bar{P}/P_\infty$.

A few slices of the pressure field obtained with the matrix omnidirectional method, comprising of $212\times 266\times 101 = 5.6$ million grid points, are presented in figures 16(a), (c) and (d). The spatial coordinates are normalized with respect to the fin spacing W = 63.5 mm. A comparison slice at $y/W = 0.2$ showing the pressure-Poisson solution with Neumann boundary conditions is also provided in figure 16(b), showing how far off the Poisson equation solution is, even though the same reference point at $\{x/W, y/W, z/W\} = \{0,0.2,0\}$ was used in both cases, being enforced to $P/P_\infty = 1$.

Zoom In Zoom Out Reset image size

Before further discussion is provided on the pressure fields shown in figures 16(a), (c) and (d) is provided, we first demonstrate that the pressure fields obtained by the omnidirectional matrix method are a good representation of the true pressures. In figure 17(a) the pressures obtained by solving the omnidirectional matrix equation are presented at the center-plane, at the closest positions from the surface of the model (∼5 mm away from the surface), and compared with the pressures measured at the surface by a set of pressure taps connected to a Kulite KMPS-1-32-10psid pressure scanner (details in Seckin et al (2023a)). The pressures from the omnidirectional matrix method (blue line) match very well with the measured surface pressures, being spatially displaced towards the downstream direction as they are inside the lambda shock structure and miss some of the boundary layer effects. The pressures from the Poisson solver, however, are not a match with our measurements. Different boundary conditions for the pressure-Poisson solver were tested (only considering the information known a priori, i.e. the free stream pressure), none of which yielded any better results. In figure 17(b), a comparison of different residual cutoffs is performed for the dual fin case for the omnidirectional matrix method. For the residual $\epsilon_R = 10^{-4}$, 1072 iterations were required, taking 27 min to converge. The most converged case tested, at $\epsilon_R = 10^{-5}$, took 2013 iterations (40 min). The convergence behavior at the boundary is a scaling of the entire pressure field towards the converged value. The same behavior was observed for the slanted cylinder flow presented in section 4.1.

Zoom In Zoom Out Reset image size

Although the time required to converge (27 min) is fairly long for one field (the average pressure field), the mesh size is far larger than was previously explored in past studies. The pressure fields displayed in figures 16(a), (c) and (d) are data that was previously unavailable for this problem coming exclusively from experimental measurements, and highlights the main features of this flow with unprecedented detail. Features such as the oblique fin-tip shocks and reflected shocks are revealed to have a three-dimensionality due to the open top surface of this double-fin geometry. An expansion fan at the more downstream locations ($x/W\sim 1.2-1.4$) is also revealed. Close to the floor surface where the shock interacts with the boundary layer, a large lambda shock is observed at the center plane. Especially for the double-fin problem, these features are difficult to see with a shadowgraph apparatus due to the difficult optical access owing to the presence of the fins.