Modified (n+1)D Laplacian for smooth pressure reconstruction based on time-resolved velocimetry (1): analysis and numerics

We analyze a smooth pressure solver based on the ‘modified Poisson equation’: ∇2p+ξ2∂2p∂t2=f(u(t)), where p is the pressure field, u(t) is the velocity field measured by time-resolved image velocimetry, and ξ2 is a tunable parameter to control the solver’s diffusive behaviour in time. This modified Poisson equation aims at obtaining smooth pressure fields from potentially noisy image velocimetry measurements, and is a part of the current four-dimensional (4D) pressure solver (implemented in, for example, DaVis 10.2) by LaVision. This work focuses on investigating three aspects of the ‘modified Poisson equation’: smoothing effect, error propagation, and drift in time. We first provide rigorous analysis and validate that this solver can sufficiently smooth the computed pressure field by setting a large enough ξ2. However, a large value of ξ2 may cause large errors in the reconstructed pressure fields. Then we introduce an upper bound on the error in the reconstructed pressure fields to quantify the error propagation dynamics. Finally, we discuss the potential drift due to the partitioning in time, which is an optional strategy used in LaVision’s current 4D pressure solver to reduce computational costs. Our analysis and validation not only show that careful choice of the parameters (e.g. ξ2) is needed for smooth and accurate pressure field reconstruction but provide theoretical guidelines for parameter tuning when similar pressure solvers are used for time-resolved image velocimetry data.


Introduction
Velocity field measured by Time-Resolved Particle Image Velocimetry (TR-PIV) enables time-resolved pressure field and force reconstruction (Beresh 2021).It is intuitive to compute the pressure field frame by frame for sequential pressure fields.However, it is also expected that this practice may lead to noisy pressure reconstruction (Charonko et al 2010, De Kat and Van Oudheusden 2012, Panciroli and Porfiri 2013).To be more specific, the noise in the pressure field usually resemble high-frequency temporal oscillations in the reconstructed field caused by the non-correlated noise in the measurement at different frames, and this oscillation is not necessarily physical.For example, figure 1 shows a typical pressure signal sampled at the pressure field reconstructed from tomo-PIV compared against the measurement from a pressure transducer in time (a) and (b) and frequency (c) domain (adapted from De Kat and Van Oudheusden (2012)).The pressure computed from PIV has high-frequency noise (figures 1(a) and (b)) that leads to overestimated power spectral density for high frequencies (highlighted by the blue shade in figure 1(c)).This phenomenon can be observed in several studies.Examples include almost all the studies surveyed by Van Oudheusden (2013) reporting power spectral density of PIV-based pressure field, and it is not limited to pressure field reconstruction based on pressure Poisson equation.This non-physical high-frequency oscillation could be unfavourable, for example, in applications where the high-frequency fluctuations may be important and aeroacoustics studies (Ghaemi et al 2012, Nickels et al 2017, Zhang et al 2018).
To overcome this issue, there are in general two strategies: either suppressing the noise in the 'source' (i.e. the experimental measurements of the velocity field, pressure gradient, or the data of the pressure Poisson equation) or implementing some smoothing mechanism, directly or indirectly, to de-noise the 'end-product', which is the pressure field as a time series.One solution in the latter category is to introduce some diffusion in time to suppress the temporal oscillation in the reconstructed pressure field.LaVision's current pressure solver (implemented starting from DaVis 10.0 released in June, 2018) employs this idea by introducing the 'modified Laplacian' (Jeon et al 2017a(Jeon et al , 2017b)).The pressure field is reconstructed by solving a modified pressure Poisson equation in the form of where ∇2 is the 'modified Laplace operator' with a term of weighted second order time derivative (ξ 2 ∂ 2 p ∂t 2 ) added to the canonical Laplacian.For example, ∇2 : for two-dimensional (2D) in space, where ξ 2 is a user-defined parameter to be tuned to adjust the performance of the pressure solver.On the right hand side of the (1), u is the velocity field, ν is the kinematic viscosity, and ρ is the density of the fluid.
The more recent version of this pressure solver in DaVis (as of March, 26th, 2023) is numerically implemented on convective frames for ξ 2 ∂ 2 ∂t 2 : where ξ 2 ∂ 2 p/∂t 2 c is the 'convective pressure', which can be estimated by the velocity field measured from velocimetry.In this sense, the pressure field is considered a flow quantity that is transported by the moving bulk fluid.Nevertheless, both (1) and ( 2) share the fundamental similarity rooted in the modified Laplace operator.
The core characteristic of ( 1) and ( 2) lays in the 'modified Laplacian': by appending an additional term of ξ 2 ∂ 2 p ∂t 2 to the canonical Laplacian, 'diffusion' in time is introduced to the Poisson equation.This practice could provide smoother pressure reconstruction if a proper ξ 2 is selected, however, at the potential cost of i) violating the pressure Poisson equation derived from the Navier-Stokes equations and ii) accumulation of error over time.
Sakib et al (2021) used the four-dimensional (4D) solver (implemented in DaVis 10.2) to reconstruct the pressure field from the velocity field of an impinging synthetic jet measured by a time-resolved tomographic PIV system.In this study, synthetic jets with different length scales and frequencies but similar exit velocities (Re ≈ 8000, based on the peak exit velocity at the orifice and the diameter of the orifice) are used to study the effect of spatial and temporal resolution on the reconstructed pressure field.This study reported that, while the 4D solver removes the high-frequency noise, the reconstructed pressure fields sometimes suffer from drift in time.In addition, the study showed that with fixed parameters of the 4D solvers (ξ 2 = 1 and the number of space-time blocks4 same as the number of velocity fields), the non-dimensional frequency of the synthetic jet affects the slope of the temporal drift of the reconstructed pressure field.The discoveries reported in this experimental work inspired the current analysis.
Later, in Sakib et al (2023), a parametric study was performed based on the same impinging synthetic jet data.Various combinations of the parameter ξ 2 , the number of space-time blocks (m), and the oscillation frequencies of the synthetic jet were chosen to conduct a series of tests.This study showed that increasing the value of ξ 2 results in a smoother reconstruction of the pressure field in time.However, the slope of the temporal drift increases with increasing ξ 2 .This study also demonstrates that a larger space-time block size lowers the slope of the temporal drift and it is possible to choose a proper combination of ξ 2 and space-time block size that minimizes the temporal drift and high-frequency temporal Figure 1.A typical sequential pressure reconstruction from PIV (red lines) compared against measurement from a pressure transducer (black lines), adapted from De Kat and Van Oudheusden (2012).(a) pressure signal time series and (b) 0.1 s subset of the full series.(c) Power spectral density of the pressure signal.The blue shade marks the overestimated high-frequency region.Reproduced with permission from De Kat and Van Oudheusden (2012).
noise in the reconstructed pressure field.Although these studies demonstrate that parameter tuning may play a significant role in PIV-based pressure reconstruction using the modified Poisson equation based 4D solver, they do not offer analytical insights into the underlying mechanism of the smoothing effect or the propagation and accumulation of errors (i.e. the temporal drift).
Other early users of this 4D pressure solver include, for example, Michaelis and Wieneke (2019), Saiz et al (2022).Based on the Lagrangian particle tracking and PIV techniques, an experimental methodology for diagnosing the aeroelastic forces was proposed in Saiz et al (2022).The experimental data derived from PIV were processed with the LaVision Davis 10.0.5 software, and the pressure fields in the small regions outside the shear layer of a non-steady wake of an oscillating flexible thin plate trailing a circular cylinder were reconstructed using a Poisson solver.The results showed the feasibility of applying an optical diagnostic approach to evaluating the aeroelastic forces and estimating the residual of the force equilibrium equation.Michaelis and Wieneke (2019) presented a comparative experimental assessment of velocity, vorticity, acceleration, and pressure calculation by using two advanced flow measurement techniques: time-resolved and multi-pulse Shake-the-Box (STB) and Tomographic Particle Image Velocimetry (Tomo-PIV).In time-resolved Tomo-PIV, pseudo Lagrangian tracking is used to obtain the material acceleration and pressure from the 3D3C velocity field.For the multi-pulse STB, the particle tracks derived from STB were processed with both binning and fine-scale reconstruction to compute the velocity, acceleration, and pressure.In terms of pressure calculation, a Poisson solver was applied to reconstruct pressure fields from velocity fields when fine-scale reconstruction was not used.The study evaluates the accuracy and resolution of each technique and highlights the importance of selecting the appropriate flow measurement technique depending on the specific application requirements.However, these studies did not report the detailed parameter settings for the pressure solvers and how the parameters would affect the error propagation in the reconstruction.
The parameter selection of the 4D pressure solver (e.g.how to choose ξ 2 and the number or the size of the space-time blocks when partition in time is practiced) remains unclear until the fundamental property of the modified Laplacian is known.Thus, careful analysis of the modified pressure Poisson equation, to be more specific, the characteristics of the modified Laplacian in the context of PIV-based pressure field reconstruction is needed.In this work, we analyze the error propagation dynamics of this pressure solver and provide insight into parameter tuning that can mitigate the temporal noise and error accumulation.
The paper is structured as follows.Section 2 analyzes the modified pressure Poisson equation from three perspectives (i.e.smoothing, error estimates, and drift due to partitioning in time).Section 3 presents the validation and verification tests.Section 5 concludes the work.In section 4, we comment on the possible optimal choice for the parameter ξ 2 .Throughout this paper, we provide some remarks on the practice for the experimentalists.

On the smoothing effect
We first provide analytical intuition for the smoothing effect of the modified pressure Poisson equation from two different perspectives (i.e.functional analysis and Fourier analysis), with each providing different insights.At the end of this section, we also point out some relevant references from the image processing community, which also leverage the smoothing effect of a Poisson solver.
A functional interpretation of the modified Poisson equation can provide some intuition for the role of ξ 2 in noise reduction.Solving a Poisson's equation, ∇ 2 p = f, is equivalent to minimizing the functional of which requires suppressing the gradients of p across the domain Ω in all directions (Evans 2010).Similarly, by adding ξ 2 ∂p 2 ∂t 2 to the canonical Laplacian and solving ∇2 p = f is equivalent to minimizing where the integral in the additional temporal dimension Θ is involved.In turn, the computed pressure field is smoothed out due to the suppressed gradient in the time component.The factor ξ 2 adjusts the relative weight of time derivative (∂p/∂t) over the gradients in space (∇p).For example, a sufficiently large ξ 2 can make ||ξ ∂p ∂t || dominate Ȇ, and minimizing Ȇ or solving the corresponding modified Poisson equation is a matter of reducing the time derivative of p alone.The resulting solution is expected to be smooth.Heuristically, we expect that larger ξ 2 would lead to smoother p in time.On the other hand, any potentially high temporal fluctuation in the solution (i.e. a high value of ||∂p/∂t||) would lead to large || ∇p|| 2 even if ξ 2 is small.This implies that a small ξ 2 could be effective in terms of smoothing p in time.Note, both the signal and noise in the solution p are smoothed in the same way.
The other way to explain the smoothing effect is to consider solving a Poisson equation as running a low-pass filter on the data f (De Kat andVan Oudheusden 2012, Faiella et al 2021).The solution of the Poisson equation (p) is almost always smoother than the original data (f) in all dimensions due to the property of investing a Laplace operator.This statement also holds for the modified Poisson solver and is independent of the specific numerical methods for the solver (Faiella et al 2021).We now use a simple example to briefly elaborate on this point.Consider a one-dimensional (1D) time-varying flow field u(x, t) on x ∈ (0, X), and we want to use the modified Poisson solver (i.e.(1)) to compute the pressure field of the flow over time t ∈ (0, T) with some boundary conditions.This task requires a two-dimensional (2D) solver (one dimension for space and one for time).Assuming the solution of (1) takes the form of where ϕ j (x) and c k (t) are eigenfunctions with corresponding eigenvalues λ j,k , and A j,k are the coefficients to be determined (j and k are integers).Plugging ( 5) into (1) and expanding f with the same basis of ϕ j (x) and c k (t) lead to where a j,k are the coefficients.The solution of (1) can be found by comparing the terms in (6) and A j,k are evaluated The coefficients A j,k can be interpreted as an amplification ratio for each frequency component of f to the same frequency component in p (see Faiella et al (2021) for more detailed illustration and examples).For any physical flow field, the values of a j,k are finite, and the high-frequency components decay.In other words, a j,k decreases with larger j and/or k in general.In addition, λ j,k grow fast with increasing j and k.For example, sufficiently large j and k leads to λ j,k ∼ −( j 2 X 2 + ξ 2 k 2 T 2 ), and A j,k diminishes for higher frequency components, both in time and space.Thus, solving the modified Poisson equation behaves like passing the data f (signal and noise) through a low-pass filter, which features the smoothing effect.One may also notice that larger ξ 2 contributes to larger λ j,k and is more effective in suppressing the 'propagation' of high temporal frequency components of c k (t) from f to p. Again, smoother p in time is expected for larger ξ 2 .The above analysis holds for problems in higher dimensions as well.We will see some demonstration of how the choice of ξ 2 affects the smoothing effect in section 3.
Remark (Tuning ξ 2 for broadband signals).For broadband signals with significant high-frequency content, if which is also of interest, one has to apply less regularization (i.e.smaller ξ 2 ) to preserve the high-frequency components.However, when trying to conserve the high-frequency signals with small ξ 2 , high-frequency noise would be retained in the same way.This renders an intrinsically difficult challenge to improve the signal-to-noise ratio over a wide spectrum, from both the experimental and pressure solver design perspectives.If the low-frequency components are of more interest, aggressive regularization with large ξ 2 is possibly appropriate.
In fact, Laplacian is widely used in the fields of imaging processing and computer graphics to sharpen images, which accentuates transitions of an image and, in turn, increases fluctuation over the domain (i.e. the image in this context).Solving a Poisson equation, which is equivalent to applying an inverse of Laplacian to the data, provides smoothed images.This smoothing effect is leveraged in many classic tasks such as blending (Pérez et al 2003), seamless image stitching (Levin et al 2004), smooth inpainting (Bertalmio et al 2001), and temporal smoothing for videos (Bonneel et al 2015).

On the error propagation
In this section, we analyze the error propagation dynamics of the modified Poisson solver.Consider that we are reconstructing a pressure field on a domain Ω in space, over a period of time in Θ, using the method described in (1).Due to the error in the data field (ϵ f ), the pressure field computed by solving the modified Poisson equation (p) is also contaminated: where ϵ p is the error in the computed pressure field.The true value of the pressure field is governed by the canonical pressure Poisson equation when the data f is not corrupted: Comparing ( 7) with (8) leads to which is a modified Poisson equation with respect to ϵ p , and εf is the total error in the data.Multiplying both sides of (9) by ϵ p and integrating by parts leads to ˆ∂Θ ˆ∂Ω If we ignore the error contribution from the boundary, the first term in the right hand side of (10) vanishes5 .We can estimate the error level in the reconstructed pressure field: error in data where C is the Poincare constant6 , which is associated with the smallest eigenvalue of the boundary value problem of (9) on a specific domain.Factors such as the dimension (i.e.1D, 2D, 3D, or 4D), shape and size, and boundary condition setup of the domain couple determine the value of C, which in turn dictates the error in the reconstructed pressure (Pan et al 2016, Nie et al 2022).The user-defined parameter ξ 2 also presents in C, in addition to serving as the 'weight' of the additional error source of in (11b).Thus, the value of ξ 2 could significantly affect the error propagation.It is also worth noting that the additional error source on the right hand side of (11b) is related to the true value of the temporal fluctuation of the pressure (measured by second order time derivative, to be more specific).We will illustrate the rich error propagation dynamics of this modified pressure Poisson solver using some examples in section 3.

Remark (The cost of smoothing).
Larger ξ 2 provides more aggressive smoothing in p, but at the cost of introducing extra error, which is associated with the power of the pressure fluctuations ||∂ 2 p/∂t 2 || L 2 (Ω+Θ) .This error is rooted in the (slight) violation of the pressure Poisson equation and is not physical.Thus, ξ 2 should be as small as possible, given that the reconstructed pressure is smooth enough.
We can also use the diffusion equation to gain some intuition about the characteristics of the modified Poisson equation.
Besides the application of the Poisson equation in the context of the pressure field, a Poisson equation also governs the steady state of a diffusion process.For example, (9) can be thought of as the steady state of the diffusion of ϵ p driven by the source of εf .Interestingly, one should not be confused with the modified Poisson equation ( 9) with the canonical diffusion equation.A diffusion equation involves the first-order derivative of time.For example, ∂q ∂t = ∇ 2 q − s is a general diffusion equation with a unity diffusion coefficient in all spatial dimensions, where −s is the source, q is the dependent variable.At steady state, the corresponding Poisson equation is ∇ 2 q = s for t → ∞.The 'arrow of time' applies here for t, while the spatial coordinates feature the symmetry in space.However, the modified Poisson equation, ∇ 2 q + ξ 2 ∂ 2 q ∂t 2 = s, augments the spatial dimensions with one more dimension for diffusion in time.Here, the time dimension has no fundamental differences from spatial ones by the mathematical expression.In other words, the added time dimension can be considered a new 'spatial' dimension, and the time t in this setup features symmetry.
The modified Poisson equation can be considered a model of the steady state of an anisotropic diffusion process in a heterogeneous media, where the diffusion coefficient in t is ξ 2 while the ones for the coordinates in space are all unity.Small ξ 2 means that q diffuses slowly in the directions of t (forward and backward in time), and q at adjacent moments can be distinct.If q represents ϵ p , as the case of (9), low temporal diffusion due to small ξ 2 implies that high-frequency fluctuations (noise or signal) in time is possible.If ξ 2 is relatively large, high diffusion in the directions of time would spread q to (immediately) subsequent and previous time instants, and thus q cannot drastically vary in time.In turn, smoothing in time is achieved by the 4D pressure solver if q stands for ϵ p .However, smearing of the true variable is also introduced when q represents p, as in (1).Tuning the value of ξ 2 allows a tradeoff between these two effects.

On the partitioning in time and drift
Solving a (n + 1)D pressure Poisson equation can be computationally expensive (more arithmetic operations and higher memory demand than the corresponding nD solver) especially when the mesh on t is large due to time-resolved measurements.One strategy to overcome this challenge is to partition the domain into smaller subdomains to localize the computation (see the partition of unity method used for example in Larsson et al (2017)).LaVisoin's 4D solver implemented similar practice (Jeon 2023), and an illustration using a (1 + 1)D flow on x ∈ [0, X] over t ∈ [0, T] is shown in figure 2(a).t is defined as the ratio of the time and the time span T of the measured data.The entire domain of X × T is divided into m spacetime blocks, and each space-time block has a duration of T/m.At t = 0 the pressure field is computed from (n = 1)D data and used as the Dirichlet boundary condition (indicated by the purple solid line at the left edge of the first space-time block in figure 2(b)).For an arbitrary time, Dirichlet boundary conditions are prescribed at the edges of the spatial domain (marked by the pink solid lines in figure 2).At the end of i-th space-time block, a Neumann boundary condition can be prescribed by leveraging the convection of pressure fields (Jeon 2023).This prescribed Neumann boundary condition is indicated by the green dashed lines in figure 2. After the reconstruction of the pressure field in the i-th space-time block, the pressure field at the end of this space-time block (i.e. the right edge of the i-th space-time block) is known.This computed pressure field is then adopted by the next space-time block as a Dirichlet boundary condition (e.g. the green solid line at the left edge of (i + 1)-th space-time block for x).
In this section, we provide a preliminary analysis based on the (1 + 1)D problem to 'qualitatively' demonstrate the drifting caused by the partitioning in time, and visualize the analysis with some examples shown below in figure 3. Let us consider the error propagation problem of ( 9) in (1 + 1)D.If we ignore the error on the boundaries for the computation in the i-th space-time block (denoted as (Ω + Θ) (i) ), ( 9) becomes where εf = ϵ f + ξ 2 ∂ 2 p ∂t 2 is the total error in the data.Here, we introduce a re-scaling of time t: τ = t/ξ is a normalized time scale.This allows us to transfer the modified Poisson equation to the familiar canonical Poisson equation and understand how the value of ξ affects the diffusion of the error.
If εf is expended as to match the boundary conditions and indicates the end of previous space-time block.The solution to (12) corresponding to the first eigenvalue 7 is 7 The first eigenvalue and the corresponding solution is often important if not dominant in a real-world problem. where is the first eigenvalue of the problem of (12).Thus, ϵ p 1,1 takes the form of λ1,1 sin(π x X ) ̸ = 0 at the right boundary of (Ω + Θ) (i) , which is end of this space-time block (t = iT/m).
In turn, at the end of the i-th space-time block, ϵ p is not necessarily zero (denoted as ϵ (i) p ), and this error will be adopted by the next space-time block as the error on the Dirichlet boundary condition.This observation provides insights into the error accumulation over time-or we may call it 'drift'-in this space-time block.We want to emphasize that the expansion terms of εf have the opposite signs to the corresponding components of ϵ p , and the signs of a j,k determine the direction of the drifting.For example, if εf is dominated by some bias error, say εf = −1, the coefficient of the dominant component is a 1,1 < 0 and ϵ p 1,1 > 0, meaning that the reconstructed pressure is expected to drift up.In fact, depending on how the data in (1) or ( 2) is formulated, negative definite bias may occur in εf (see Nie et al (2022) for more details).If this is the case, it is almost guaranteed that ϵ (i) p > 0 and upward drift would occur in the computed pressure field.
Moving to the next space-time block, the error propagation in the (i + 1)-th space-time block can be modeled similarly, however, by a non-homogeneous problem: The solution to this problem can be found by superposing a homogeneous problem similar to ( 12) and a nonhomogenous Laplace equation: p takes the form of ϵ p 1,1 , the solution to (15) is where A = exp(− π T ξ mX ) and B = exp( π T ξ mX ).This solution can be thought of as the error on the 'left' boundary (t ).The full solution to ( 14) is a superposition of ( 16) and If the domain of (Ω + Θ) (i+1) is long enough in the time dimension (physically requires small ϵ and/or small m for a given T), the error on the boundary at t = i T m , which is adopted from the i-th space-time block, is dissipated and vanishes towards the end of the (i + 1)-th space-time block.If the domain of (Ω + Θ) (i+1) is not long enough and does not allow sufficient decay of ϵ p -by the end of the (i + 1)-th spacetime block, the residual error at t = (i + 1) T m from the Laplace equation overlays the accumulated error from the Poisson equation, which is the same as ϵ (i) p .The total error in the pressure field at t = (i + 1) T m , denoted as ϵ , for the (i + 1)th space-time block is thus expected to be higher than ϵ The same monotonic error growth would occur throughout the entire domain of (Ω + Θ), from the first space-time block to the last one.Thus, ϵ (1) and drifting is expected if there exists a dominating negative definite component in εf .
We now use various combinations of ξ 2 and the number of space-time blocks m to visualize the above analysis.Figures 3(a) and (c) shows the error profile according to ( 12) and ( 14) with ξ 2 = 1 and ϵ f = −1.It can be observed that the error in the (i + 1)-th space-time block is higher than the error in the i-th space-time block and the reconstructed pressure is expected to drift up.Figures 3(a)-(c) illustrates the superposing of the homogeneous problem ( 12) and the nonhomogenous Laplace equation (15) for the (i + 1)-th space-time block.The reconstructed pressure of the (i + 1)-th space-time block drifts up because the domain of (Ω + Θ) (i+1) is not long enough in the time dimension for error dissipation.Compared figure 3(d) with figures 3(a) and (c), we can see that the number of space-time blocks affects the drift.With a larger number of space-time blocks m, it takes a longer time for the error to grow monotonically to the maximum value as shown in figure 3(c).
Additionally, comparing figures 3(e) and (f) with figures 3(a) and (c) shows that the value of ξ 2 also has Figure 3. Error in the i-th space-time block (a), with relatively large weighting factor ξ 2 = 1.Superposition of the numerical solution of ( 12) and ( 15) yields to the solution of ( 14), i.e. panels (a) + (b) = (c), where (c) is also the error in the (i + 1)-th space-time block.Error in the reconstructed pressure field with a larger number of space-time blocks, m = 6 in the time span of 0 ⩽ T ⩽ 3, (d).The error in i-th and (i + 1)-th space-time blocks with a smaller weighting factor ξ 2 = 0.1, (e) and (f).
an effect on the drift, given a fixed number of the spacetime blocks (m = 2 in this case).In figures 3(a) and (c), the error reaches the maximum value at the beginning of the (i + 1)-th space-time block.However, figures 3(e) and (f) shows that with the decrease of ξ 2 , the error accumulates faster and approaches the maximum early in the i-th space-time block.
Invoking the diffusion analogy in section 2.2 to (15) may provide an alternative explanation for the drift.The error accumulated by the right end of i-th space-time block takes 'time' in the next space-time block to dissipate (by the boundaries).If ξ 2 is large, high diffusivity will effectively transfer the error on the left end of the (i + 1)-th space-time block to the right boundary.While small ξ 2 can isolate the error propagation in the directions of time.For a given T, large m means that each space-time block is short in the temporal dimension, and the 'residual' error on the left boundary of the (i + 1)-th spacetime block can be easily conducted to the other end in time, which will be picked up the next time domain.Thus, using small ξ 2 and fewer space-time blocks may eliminate the drifting.Recalling that the above analysis is based on the assumption that there potentially exists a 'direct current' bias in εf due to a formulation of the data commonly used in many works (see for example Gurka et al (1999), De Kat and Van Oudheusden (2012)).Using an alternative formulation of the data (Nie et al 2022) may also be helpful in eliminating the definite drift.
Remark (Beware of the drift).Diffusion in time introduced by the modified Poisson solver may cause drift in the reconstructed pressure field.The drifting behavior can be complex and tuning m and ξ 2 may be able to control the drift depending on the inherent nature of the flow, error, as well as the setup of the domain.Without a priori data and analysis, it is necessary to be attentive to the reconstruction.For example, comparing the reconstructed pressure against the measurement by a pressure transducer is a simple and effective way to monitor the potential drift.

Validation on (n + 1)D examples
In this section, we use (n + 1) dimensional, denoted as (n + 1)D modified Poisson solvers, to illustrate the core features of the original 4D pressure solver, where the solver has n-dimensions in space with one extra dimension in time.(Charonko et al 2010) in a channel with a length scale of X and time scale of T, driven by a dimensionless time-varying pressure gradient of ∂p/∂x = 1 + γ sin(π kt/T), where γ is the ratio between the strength of the steady and oscillating pressure gradients, t is the dimensionless time, x is the dimensionless coordinate along the channel and k is the wave number.This leads to an one dimensional (1D) pulsatile flow on x ∈ [0, X], and the Laplacian of the pressure field is The corresponding boundary conditions can be: p(x = X, t) = 1 + γ sin π kt T X and p(x = 0, t) = 0. Integrating the pressure gradient along x directly gives the pressure field (see figure 4(a) for an illustration): To solve ( 19), two other boundary conditions for t = 0 and t = T are needed.The most straightforward idea is perhaps to directly adopt the pressure field computed from the canonical Poisson equation (i.e. ( 17)).Thus p(x, t = 0) = x and p(x, t = T) = [1 + γ sin(π k)]x are known and can be used as Dirichlet boundary conditions.For such a setup, if we assume the error on boundaries is negligible and only consider how the error in f affects the computed pressure field, the problem is manageable and rich enough to expose the key properties of (1).
For the setup of the canonical Poisson equation ( 17), the smallest eigenvalue is λ 1 = π 2 X 2 , and the corresponding Poincare constant is The upper bound of the error in the pressure field is: for each frame and the time average of the error bound takes the same form.For the setup of the modified Poisson equation ( 19) the space and time-averaged error in the computed pressure field takes the form of (11b), and the Poincare Comparing ( 20) and ( 21), we expect distinct error propagation dynamics from data to the reconstructed pressure field (ϵ f → ϵ p ) for the canonical nD pressure solver and (n + 1)D solver.We next illustrate the key features (i.e.smoothing effect and error propagation) of these two different solvers using an example with specific parameters.
We first discuss the impact of ξ 2 on the error in the pressure field.We used the centerline of the viscous pulsatile flow From figures 4(g)-(i), we can see that with increasing ξ 2 (from 0.001, 0.1 to 0.1), the profile of the error in the reconstructed pressure field is smoother, however, the amplitude of the error in general increases, resulting higher power of the error (∥ϵ p ∥ L 2 (Ω+Θ) = 0.0898, 0.0902, and, 0.1717).In this particular case, we can see that a proper choice of ξ 2 can reduce the error in the pressure field (e.g.ξ 2 = 0.001 leads to smaller ∥ϵ p ∥ L 2 (Ω+Θ) than the results from 1D solver and (1+1)D solver for ξ 2 = 0.01 and 0.1), in addition to the smoothing effect.
This phenomenon may be more clearly observed by sampling the pressure at a certain location in the spatial domain over time.5 shows that with larger ξ 2 , the reconstructed pressure is increasingly smoother, but a higher error may be expected.This example demonstrates that the value of ξ 2 should be selected properly when using a (1 + 1)D pressure solver to strike a balance between smoothing and  error reduction.This claim also holds for higher dimensional problems.
Our analysis (i.e. ( 11b) and ( 21)) indicate that in addition to ξ 2 , the inherent frequency of the pressure fluctuations of the flow (measured by ∥∂ 2 p/∂t 2 ∥ L 2 (Ω+Θ) ) also impacts the error propagation.For example, in the pulsatile flow we used here for demonstration, ∂ 2 p/∂t 2 ∼ k 2 and the value of wave number k affects the power of ∥∂ 2 p/∂t 2 ∥ L 2 (Ω+Θ) .We chose k = 3, 9, and 15 to illustrate the impact of k on the error propagation.Parameters other than k are adopted from the value we used in figure 4(a).The exact pressure field with different k are shown in figures 6(a)-(c).The computed pressure field and the corresponding error in the pressure field are illustrated in figures 6(d)-(i).It is observed from figures 6(g)-(i) that when the fluctuation frequency of the pressure field is higher, using (n + 1)D pressure solver may produce a higher error (e.g.∥ϵ p ∥ L 2 (Ω+Θ) = 0.0895, 0.0902, 0.1077 for k = 3, 9, 15, respectively).The phenomenon may be more clearly revealed by the reconstructed pressure field and corresponding error sampled at x = 0.5 (shown in figure 7).Inequality (21) also indicates that the upper bound of the error is significantly affected by the length scale (X) and the time scale (T) of the domain.This is similar to the impact of the geometry of the domain (i.e.size and aspect ratio) on the error propagation as analyzed in Pan et al (2016).Here, we chose four combinations of various X and T with other parameters fixed (i.e.ξ 2 = 0.01, ϵ f = −1 + e, e ∼ N (0, 1) and k/T = 3) to test how the error propagation dynamics is affected when the geometry of the domain is changed.Figure 8 illustrates the pressure fields for different combinations of X = 1 or 3 and T = 1 or 3.The relative error (ϵ p /P max ) in the pressure field reconstructed from the (1 + 1)D solver is shown in figure 9, where P max is the maximum pressure in the domain.
In figure 9, it can be observed that among the four cases, the error in the pressure in figure 9(d) is the highest, followed by the case in figures 9(b) and (c), and the error of the case in figure 9(a) is the lowest.The quantitative measure of the error for each case, ∥ϵ p ∥ L 2 (Ω) , confirms the above observation (see the caption in figure 9) and is consistent with the prediction by (21).We next move on to analyzing the effect of parameter ξ 2 on the asymptotic behaviour of error propagation.When ξ 2 → 0, (21) reduces to (20), which is equivalent to a 1D solver as expected, and smoothing effect vanishes.When ξ 2 is sufficiently large (ξ → ∞), (21) transforms to Interestingly, this indicates that the error in the reconstructed pressure field is dictated by the additional error term, ∂ 2 p/∂t 2 , introduced by using an 'over-diffused' (n + 1)D solver.The Poincare constant (C = T 2 /π 2 ) in ( 22) is a sole function of T, which means that the total error is dictated by the length of the temporal dimension (T) alone and the spatial length scale of the domain (X) does not matter.Despite that large ξ 2 suppresses the fluctuation in the reconstructed pressure ('signal' or noise), the solution may not be physical anymore.Therefore, a trade-off between smoothness and accuracy is needed by choosing a suitable value for ξ 2 .However, this may not be a trivial task even if the error on the boundary is ignored, as the geometry (e.g. the relative scale of X and T) and type of boundary conditions of the domain, the relative power of the two error sources (ϵ f and ∂ 2 p/∂t 2 ) coupled with the parameter ξ 2 determine the error propagation from the data to the reconstructed pressure field.We provide some preliminary remarks about selecting optimal ξ 2 in section 4.
Normalizing (21) by the power of the error in the data leads to the upper bound for the ratio of the error in the reconstructed pressure field to the error in the data: where is the ratio of the flow fluctuation to the error level in the data.Ar * (ξ) = ∥ϵ p ∥ L 2 (Ω+Θ) /∥ϵ f ∥ L 2 (Ω+Θ) can be interpreted as an error amplification ratio.
We carried out numerical experiments by solving the modified pressure Poisson equation with artificial error introduced (i.e. ( 7)) based on the aforementioned 1D pulsatile flow as shown in (18).For various combinations of X, T, and F, we calculated Ar * for different ξ 2 from the numerical experiments and compared them with the analytical bounds.In figure 10(a), we show that the values of Ar * obtained from the numerical experiments (represented by square symbols in various colors) well agree with the theoretical prediction given in ( 23), denoted by the solid curves with corresponding colors.
The above validation shown in figure 10 is performed on the basis of flows with k = 1.However, The error bound calculated by the Poincare constant based on the smallest eigenvalue may overestimate the real error, especially when high-frequency components are involved in the data.For example, in the context of the pulsatile flow we used here for validation, k can vary.An improved error estimate can be achieved by taking k into account for the eigenvalues and the error in the reconstructed pressure field can be estimated as We conducted additional numerical experiments using a flow that pulsates at a higher frequency (e.g.k = 3) to validate (24).Figure 11 provides a comparison between the values of Ar * obtained from these experiments (represented by square symbols), the error bound (depicted by solid curves based on (23)), and the error estimate considering the higher frequency in the flow (represented by dashed lines in corresponding colors based on ( 24)).The data points from the numerical tests are observed to lie below the upper bound defined by (23).However, when ξ 2 is large, (23) tends to overestimate the error.On the other hand, the improved error estimate indicated by the dashed lines aligns better with the trends exhibited by the data points from numerical experiments, despite not being strict bounds for the error as in ( 23).The error estimate derived from ( 24) is valuable in predicting the trends of Ar * as a function of ξ 2 for different frequency components in the data.As demonstrated in figure 11, depending on various factors such as the experimental setup (e.g. the selection of X and T) as well as the quality of the experiments and the characteristics of the flow field (e.g.F and k), Ar * may either increase or decrease with varying ξ 2 .Consequently, it becomes apparent once again that a careful choice of ξ 2 that can suppress the error propagation is needed.This claim holds for other boundary condition setups (e.g.see appendix A for a (1 + 1)D example with mixed boundary conditions) and higher dimensions (see, for example, a (2 + 1)D demonstration in appendix B).
Figures 10 and 11 imply that k, associated with the temporal frequency of the flow, is also important for the error propagation of the modified Poisson equation.This concurs with the experimental observation in Sakib et al (2023).To demonstrate the influence of k on Ar * , invoking (18) as an example, we can explicitly find the norm of the second time derivative of pressure field 6 , which can be substituted into (24).Inequality (24) then becomes: and the error amplification ratio turns into: Note, ( 26) indicates that the error amplification ratio depends on the temporal frequency of the flow.Assuming a unity error level in the data (i.e.∥ϵ f ∥ L 2 (Ω+Θ) = 1) Figure 12.Ar * as function of k for various combinations of X, T, and ξ 2 for the Dirichlet case.The solid lines (calculated with the smallest eigenvalue) and dashed lines (calculated with varying k) represent the theoretical upper bounds based on inequality (26).The data points illustrate the simulation results when k is changed.
for k → 0, and for high frequency flows as k → ∞.
We conducted numerical experiments using the (1 + 1)D solver and based on the pulsatile flow with different combinations of X, T, and ξ 2 , and varied k to validate (26) and illustrate the impact of flow frequency on the error propagation.Figure 12 compares the results from the numerical experiments and the estimates based on (26).The solid lines represent the upper bounds by the optimal Poincare constant which is associated with the smallest eigenvalue (i.e.Ar * (k) = , and the wave number k in the data varies from 0 to 30).The dash lines are derived from the error estimation based on (26) counting for varying frequencies of the flow.It can be observed that the error estimates (dash lines) agree fairly well with the trends shown by the numerical experiments (square symbols with corresponding colors), and they are enclosed by the upper bounds (solid lines) based on the optimal Poincare constant.However, the upper bounds overestimate Ar * (k) for high-frequency flows (i.e.k is large).When k increases, Ar * (k) asymptotically plateaus at constants that depend on X (recall ( 28)).In figure 12, we can observe that the test groups with the same value of X yield Ar * (k) converge to the same constant as predicted by (28).

Partitioning in time and drift.
To validate the analysis for the drifting effect discussed in section 2.4, we solved the (1 + 1)D modified Poisson equation, with various combinations of m, ξ 2 and εf , to recover the pressure field of the aforementioned flow as shown in figure 4(a).Figure 13 shows the reconstructed pressure fields and the error in pressure fields by solving the (1 + 1)D Poisson equation with partitioning in time (with εf = −1 being fixed while varying m and ξ 2 ).Compared figure 13(b) with figure 13(d), it can be observed that the number of space-time blocks (m) affects the error accumulation (e.g. the maximum error and the time required for the error to drift to the maximum value).In terms of the value of ξ 2 , it can be seen by comparing figures 13(b) and (f) or (d) and (h) that when m is fixed, with the larger value of ξ 2 , the maximum error becomes lower and the time required for the error to drift to the maximum becomes longer due to the effective diffusion with higher value of ξ 2 (seen in section 2.4).
Figure 14 shows the error in the reconstructed pressure field at x = 0.5.It can be observed from figure 14(a) that with fixed values of ξ 2 and εf (i.e.ξ 2 = 0.01, εf = −1), long term and slow drift may occur for a large number of partitioning spacetime blocks while using a smaller number of space-time blocks lead to a fast drift.Over time, the drift stops, and the error oscillates near a constant bias.Moreover, figure 14(b) shows that ξ 2 also affects the drift (m = 100, εf = −1 are fixed in this test).Slower drift occurs with a larger value of ξ 2 , and over time, the drift vanishes.This observation resembles the test shown in figure 14(a).This similar effect of m and ξ on the temporal drift perhaps can be explained by noting the normalization in (12).The time scale of each space-time block is τ = T ξ m , where ξ and m are interchangeable.Last, figure 14(c) demonstrates that the direction of the drift can be affected by the profile of the error in the data ϵ f (tests based on fixed values of m = 400, and ξ 2 = 0.1 while varying ϵ f ).By this 'minimal' demonstration (e.g. with simple errors in the data field, but without any error on the boundaries), one may already tell that the drift behavior is complex.
To further explore the drifting behavior when partitioning in time is used for (n + 1)D, we test the (1 + 1)D solver with a different boundary condition setup on a domain of x ∈ [0, X] over t ∈ [0, T] as shown in figure 15(a).Dirichlet boundary condition (computed or adopted) is only applied on the first frame of each space-time block, and the remaining boundaries of each space-time block are set as Neumann conditions.More realistic artificial error εf = 0.1 + 0.1 × ξ 2 was introduced into the data field f, where the first component εf = 0.1 represent the bias error in the data, and the second εf = 0.1ξ 2 models a simple scenario of the error caused by smoothing.This setup is inspired by ( 9) and (11b).We also introduce bias error εf = −1 to the Neumann boundaries to resemble a relatively high error    on the Neumann boundary conditions.This is often the case, as when Neumann conditions have to be applied (e.g. a wall as a boundary or an edge of the domain is in the rotational region of the flow), the error in the image velocimetry is usually high.Based on the above setup of the boundary conditions and error, we varied the combinations of m and ξ 2 to reconstruct the pressure field and test their effects on drift.
Figure 16 shows the error in the reconstructed pressure field sampled at x = 0.5.It can be observed from figure 16(a) that the drift of error becomes more significant when large m is applied (e.g.comparing the blue and red curves for m = 500 and 200 in figure 16(a)), and when m is sufficiently small, the drifting effect vanishes.In addition, figure 16(b) shows that a larger value of ξ 2 leads to more significant drift (e.g.comparing the blue and red curves for ξ 2 = 1 and 0.1 in figure 16(b)).Despite the simple test based on (1 + 1)D problem, the numerical results shown in figure 16 captures the major observations (i.e.large m and/or ξ 2 causes high drift) by Sakib et al (2023), which is based on time resolved 3D experimental data processed by the (3 + 1)D solver implemented in DaVis 10.2.Again, m and ξ show similar effects, due to the interchangeable feature indicated by the normalization in (12).
We want to emphasize that the results shown in figure 14(e.g.small m or ξ 2 causes high drift) with the tests shown in figure 16 where large m or ξ 2 leads to high drift are drastically different.This difference is caused by a 'simple' change of the boundary condition setup and the error introduced in the data.It is impossible to exhaust all setups for boundary conditions and error, etc.However, based on the heterogeneous diffusion analogy and the concert examples in this work and Sakib et al (2023), we can conclude that (i) drift may happen and its behaviour is not trivial.The drifting behaviour associated with the partitioning in time can be affected by the combined effects of boundary condition setup, the profile of the flow and error, as well as parameter tuning for ξ 2 and m; (ii) Given a data set and setup of the computation, one can tune ξ 2 and m to control the potential drift; iii) similar effects of m and ξ provide some flexibility for parameter tuning.

Remark (Prioritize velocimetry quality over parameter
tuning for ξ 2 and m).Using the (n + 1)D solver is a regularization strategy given that the velocimetry measurement is noisy.The performance of the (n + 1)D solver not only depends on ξ 2 and m, but the intrinsic nature of the flow and the error in the measurement.However, a better idea than selecting proper ξ 2 and m perhaps is to improve the velocimetry quality, which allows choosing small ξ 2 , or even avoiding using (n + 1)D solver.Thus, an experimentalist should prioritize velocimetry quality, and worry about the hyperparameter tuning for ξ 2 and m later, providing that an (n + 1)D solver would have to be used.

(2 + 1)D example
Based on the (1 + 1)D tests in section 3.1, the core features of the modified pressure Poisson solver are validated.In this section, we used (2 + 1)D pressure Poisson solver to show that the properties of the modified Poisson solvers in higher dimensions are consistent with the ones demonstrated in (1 + 1)D.Validation of the (2 + 1)D solvers based on 2D pulsatile flow used in (Charonko et al 2010) can be found in appendix B and tests based on a more complex flow are demonstrated below.
High-resolution simulations of a 2D synthetic jet impinging on a solid wall were used to test the properties of (2 + 1)D modified Poisson equation to validate our theories.An orifice with a width of 4 mm located at the bottom of a 2D rectangular domain (100 mm × 60 mm) blows a pulsing water jet at 0.5 Hz.The boundary conditions were set to be similar to the experimental setup conducted by Sakib et al (2023).A sinusoidal velocity profile was applied at the orifice as a velocity inlet, with maximum velocity U max = 0.125 m s −1 .A wall with a width of 40 mm was placed above the orifice.The free surfaces next to the impinging wall were set as zero pressure outlets.The remaining boundaries of the domain were set as walls (see figure 17(a)).As the Reynolds number of the synthetic jet based on the Stokes boundary layer thickness is small (Re δ = 100), we simulated the laminar synthetic jet at a resolution of 400 × 400 mesh in the central domain with time step of 0.005 s.After careful mesh independence check and verification, five cycles of the pressure and velocity fields were extracted from the simulation and used as the 'ground truth' to generate the synthetic data and study the behavior of the solver.
We chose the center region (a 40 mm × 40 mm window) in the flow field as the spatial domain to test the (2 + 1)D solver.Gaussian noise ϵ u and ϵ v (ϵ u,v = e, where e ∼ N (0, σ 2 ) and σ/U max = 2.4 × 10 −3 ) were introduced into the velocity field u and v, respectively.The pressure field was recovered by solving the modified pressure Poisson equation in (2 + 1)D using the contaminated velocity field: with various ξ 2 (i.e.ξ 2 = 0.001, 0.01 and 0.1).This formulation of the right hand side of ( 29) is according to a popular practice used in several studies (Gurka et al 1999, De Kat andVan Oudheusden 2012).The reconstructed pressure field p (figure 18) and corresponding error ϵ p (shown in figure 19) from solving (29) are non-dimensionalized by dividing 1 2 ρU 2 max (denoted as p * and ϵ * p , respectively).The dimensionless time t * in this section is defined as the ratio of the dimensional time to the period of the synthetic jet, and the length scales are normalized by the effective stroke length L 0 (e.g.x * = x/L 0 , where L 0 = 80 mm), which is the distance traveled by a slug of fluid during the effective forward stroke of the oscillating jet Sakib et al (2023).It can be observed that there is always bias error in the computed pressure field and the amplitude of the error in general increases with larger ξ 2 (from 0.001, 0.01 to 0.1).
Figure 20 shows the sample of the reconstructed pressure In figure 20, we can observe directional bias in the reconstructed pressure field when (29) is used.This directional bias in p * is caused by a directional error in the data field ϵ f (illustrated in section 2.4), which is from the square terms, such as (∂ϵ u /∂x) 2 and (∂ϵ v /∂x) 2 , when evaluating the data in (29).Therefore, one strategy to avoid the directional error in the  ∂y + ( ∂ṽ ∂y ) 2 = 0, and (30) only contains cross terms.As long as the noise in different velocity components is not strongly correlated or biased, the numerically evaluated data in (30) suffer less from directional bias compared to (29).This change in the formulation of the data could also change the error propagation behavior of the (modified) Poisson equation, from the velocity field to the pressure field.
We carried out a numerical test similar to the one presented in figures 18-20 by reconstructing the pressure field based on (30) with various ξ 2 (i.e.ξ 2 = 0.001, 0.01 and 0.1).Gaussian noise (ϵ u,v = e, where e ∼ N (0, σ 2 ) and σ/U max = 8 × 10 −3 ) were introduced into the velocity field as the artificial error.The reconstructed pressure field and the corresponding error from solving (30) are shown in figures 21 and 22.Compared with the figure 19, by solving (29), the error in the pressure field reconstructed by (30) (see figure 22) is lower and the bias is mitigated, despite that higher noise is introduced in the velocity field.Figure 23 shows the temporal sample of the reconstructed non-dimensional pressure field and error at the center of the domain.It can be observed that the reconstructed pressure and corresponding error are increasingly smoother  with larger ξ 2 , and the directional bias or drift is limited by using (30) for pressure field reconstruction.
Inspired by the aforementioned observation, we investigated a strategy to reduce the drift.Using (30) to eliminate the directional bias in the data could reduce the drift in time when partitioning in time is practiced.To confirm the proposed strategy, various combinations of ξ 2 and the number of space-time block m are selected to reconstruct the pressure field based on ( 29) and (30).Artificial error (ϵ u,v ∼ N (0, σ 2 ), σ/U max = 8 × 10 −3 ) was introduced into the velocity field.Figure 24 shows the reconstructed pressure sampled at the center of the domain.It can be observed from figure 24 that the number of space-time blocks m and the value of ξ 2 affect the drifting behavior of the pressure field reconstructed by (29) (represented by the solid lines), while no drifting effects present in the pressure fields reconstructed by (30) (indicated by the dash lines).Here, we want to emphasize again that the drifting behavior is complex and can be affected by combined effects of boundary condition setup, parameter settings for ξ 2 and m, flow profile, as well as the formulation of the data of the (modified) pressure Poisson equation.

Remark (Eliminate bias).
Any bias error in the data for the (n + 1)D solver could significantly corrupt the pressure reconstruction.However, random noises, in space and time, are not as destructive due to the diffusive nature of the modified Poisson equation.One should put effort in reducing the bias in the data, through high-quality velocimetry measurement and/or careful post-processing and regularization.

On the optimal choice of ξ 2
It may be possible to find an optimal value of ξ 2 so that the error in the pressure is minimized and we list several preliminary ideas here.
The first idea is inspired by the observation that the right hand side of (11b) has a minimum value by the inequality of arithmetic and geometric means: The equality of (31) holds for , and the optimal ξ 2 is This means that when the error induced by smoothing is comparable to the error in the experimental data, the optimal (or an overall balanced) reconstruction is expected.This estimation is demanding, as a priori knowledge about both the uncertainty of the error in the data (∥ϵ f ∥ L 2 (Ω+Θ) ) and the true physics of the flow field (i.e. , which is the fluctuation of the pressure in time) are needed.In addition, this ξ 2 opt would minimize the upper bound of the error in the pressure field, but not exactly the error in the pressure field.
Another attempt to find optimal ξ 2 is to directly compute the derivative of the upper bounds of Ar * (e.g. the error bound depicted by ( 23)) with respect to ξ 2 and set it to zero: It can be derived that FT 2 = X 2 is needed to minimize the global error amplification.This implies that the monotonicity of Ar * does not depend on ξ 2 but is affected by the combinations of F, T, and X.In other words, tuning ξ 2 may not be the most effective way to control the error propagation (in terms of upper bound).Instead, the boundary condition setup, design of the experiments (choice of T and X), and the control of the error (F) should be considered carefully.However, it is worth noting that the analysis above can only provide a rough guideline based on the upper bound of the error amplification ratio in the (1 + 1)D case.The actual error propagation and the optimal ξ 2 depend on the specific flow profile and setup of the solver and experiment.Additional information may be inevitable for parameter tuning.For example, the other idea demonstrated in Sakib et al (2023) leverages on the measurements in the frequency domain: ξ 2 is tuned so that the reconstructed pressure field at the location of the transducer matches the power spectral density measured by pressure transducer(s).This is perhaps a more practical solution for parameter tuning.

Conclusions and future work
We carry out theoretical analysis on the (n + 1)D modified pressure Poisson equation, inspired by the observations reported in Sakib et al (2021), which showed that the 4D pressure solver implemented in DaVis 10.2 can remove the highfrequency temporal noise, but the reconstructed pressure fields may suffer from drift in time depending on the tuning of parameters.In the current work, we analytically examined the fundamental characteristics of the modified pressure Poisson solver, including its smoothing effect, error propagation, and potential drift when partitioning in time is practiced.The analytical foundation of the smoothing effect of the modified Poisson equation is established from both functional analysis and Fourier analysis perspectives.The modified Poisson equation can be considered as adding an additional diffusion term to the canonical Poisson equation in the direction of time, and the diffusivity is controlled by a parameter ξ 2 .A smoother pressure field in time is expected for larger ξ 2 , but it may cause a higher error level of the reconstructed pressure field thus the value of ξ 2 should be selected properly.
The analysis is validated by synthetic data based on flows in (1 + 1)D and (2 + 1)D.We illustrated that the error in the reconstructed pressure field can be bounded and the error propagation of the modified Poisson equation is complex due to the intrinsic factors of the flow (e.g. the inherent frequencies of the pressure fluctuations of the flow), the setup of experiments and solver (e.g.boundary conditions setup and the geometry of the computational domain), as well as the choice of the parameter setting of the (n + 1)D solver (e.g.ξ 2 ).In addition, it is validated that the error accumulation is affected by two tuning parameters (i.e. the weighting factor ξ 2 of the temporal term and the number of space-time blocks m) when the partitioning in time is practiced.Caution must be used when choosing these parameters, and more importantly, the value of the parameters should be reported when the pressures are reconstructed using the 'modified Laplacian' based (n + 1)D pressure solvers.
The simple synthetic flows we used for validation are extendable to complex flows in higher dimensions.The (1 + 1)D and (2 + 1)D synthetic flows only involve one dominant frequency in time.This choice makes it easy to expose a fundamental property of the modified (n + 1)D Poisson solver: how the intrinsic frequency of the flow affects the error propagation.However, our analysis is not limited to a specific frequency in a flow or a flow with a dominant frequency; instead, each frequency component in a broadband signal would respect the scaling laws we developed in the analysis.Our theoretical development and core conclusions, validated in (1 + 1)D and (2 + 1)D by the synthetic data, hold in higher dimension for real experimental data (see Sakib et al (2023), where time-resolved volumetric experimental data were used to reconstruct the pressure field by the 4D pressure solver implemented in DaVis 10.2, and more comprehensive results will be presented in the part two of the collaborative series work).
The current work not only provides rigorous analysis and validation for (n + 1)D pressure Poisson solver based on modified Laplacian, but gives intuitive guidance for parameter tuning.For example, sufficiently large, but not too large, ξ 2 can effectively smooth the pressure field in time.Parameter tuning is needed to optimize the performance of the pressure solver and there is hope to establish systematic strategies for selecting a proper parameter of ξ 2 .However, this task is not trivial, and we will leave the optimization of ξ 2 and the selection of other parameters for future studies.
after normalization.
To validate the error bound in (34), we carried out numerical experiments by solving the modified pressure Poisson equation with artificial error introduced based on 2D pulsatile flow, which is similar to (18) (i.e.p(x, y, t) = 1 + γ sin π kt T x).For various combinations of X, Y, T and F, we calculated Ar * for different ξ 2 from the numerical experiments and compared them with the analytical bound.In figure 26, we show that the values of Ar * obtained from the numerical experiments (square symbols in various colors) are bounded by the theoretical prediction given in (34), denoted by the corresponding solid curves.Moreover, the values of Ar * from numerical experiments align well with the trends depicted by (34).
Taking flow frequency into account, the error estimate is Additional numerical experiments are conducted by using a flow with a higher oscillation frequency to validate (35).Similar to the results of (1 + 1)D case shown in figure 11, figure 26 shows that the data points obtained from the numerical experiments are bounded by (34), and the improved error estimate based on (35) has a better alignment with the trend of the data points.To compare the effect of spatial dimension on error propagation, we demonstrate the theoretical upper bound of Ar(ξ 2 ) * for the (1 + 1)D and (2 + 1)D cases according to ( 24) and ( 35), respectively, with various combinations of X, Y, T, F, and k (see figure 27).It can be observed that the value of the amplification ratio for the (2 + 1)D case is smaller than that of in (1 + 1)D.The (1 + 1)D case can be thought of as a special case in the limit of Y → ∞, suggesting that Dirichlet boundary condition in y direction is not available to tame the error propagation from contaminated data to the reconstructed pressure field, leading to higher error level in the computational domain.
Due to the Neumann boundary condition at t = T, we also expect higher error towards the end of the time as a Poisson equation is more sensitive to the error near the Neumann

Figure 2 .
Figure 2. The boundary condition setup on a X × T domain for a(1 + 1)D problem (a) and a partitioning in time strategy (b), which divides the whole domain into m space-time blocks.t = iT/m is the right boundary of the i-th space-time block and also the left boundary of the (i + 1)-th space-time block.

Figure 4 .
Figure 4. Pressure field and the corresponding error in the pressure fields.(a) ground truth from exact solution; (b and c) pressure and error field from 1D Poisson solver; pressure (d)-(f) and error (g)-(i) fields from (1 + 1) D pressure solver with various ξ 2 , where ξ 2 = 0.001, 0.01, and 0.1, respectively.
(γ = 1, k = 9) on a domain with a length scale of X = 1, and reconstructed the pressure field over a time span of T = 3 using the 1D and (1 + 1)D pressure Poisson solvers.Artificial error containing bias and random components (ϵ f = −1 + e, where e ∼ N (0, 1) is an independent and identically distributed point-wise Gaussian noise) were introduced into the data f to simulate the error in the experiments The results are summarized in figure 4. The exact solution of the pressure field (i.e.(18)) without contamination is demonstrated in figure 4(a), which is used as the ground truth to assess reconstruction results.Figures 4(b) and (c) show the reconstructed pressure field and the error in the pressure field (with an error level of ∥ϵ p ∥ L 2 (Ω) = 0.0919) by solving the 1D Poisson equation, respectively.We chose various values of ξ 2 (i.e.ξ 2 = 0.001, 0.01, and 0.1) for the (1 + 1)D pressure solver, and the computed pressure field as well as corresponding error in the pressure field are illustrated in figures 4(d)-(i).
Figure 5(a) shows the pressure sampled at x = 0.5 compared against the exact value (equivalent to sampling along x = 0.5 in figures 4(a), (b) and (d)-(f)).The corresponding error in the pressure is shown in figure 5(b), which is equivalent to sampling along x = 0.5 in figures 4(c) and (g)-(i).Again, figure

Figure 5 .
Figure 5. Sample of pressure (a) and the corresponding error in the pressure (b) at x = 0.5.Zoomed-in views of a portion of the samples are shown in the insets.

Figure 7 .
Figure 7. Pressure (a) and the corresponding error in the pressure (b) sampled at x = 0.5 for different k.Solid lines represent the pressure and error calculated by the (1 + 1)D solver, and dashed lines represent the ground truth of the pressure field.

Figure 10 .
Figure 10.Ar * as function of ξ 2 for combinations of X, T, and F. The solid lines represent the theoretical upper bound based on inequality (23).The data points with corresponding colors represent the computed Ar * from the numerical experiments with k = 1.

Figure 11 .
Figure 11.Ar * as function of ξ 2 for various combinations of X, T, and F. The solid lines and dashed lines represent the theoretical upper bound and estimates based on (23) and (24).The square symbols illustrate the calculated Ar * from numerical experiments with k = 3.

Figure 14 .
Figure 14.Error in the pressure ϵp sampled at x = 0.5 for various values of m-the number of space-time blocks (a) and ξ 2 (b), as well as the error in the data ϵ f (c).

Figure 15 .
Figure 15.The boundary condition setup on a X × T domain for a (1 + 1)D problem (a) and a partitioning in time strategy (b), which divides the whole domain into m space-time blocks.t = iT/m is the right boundary of the i-th space-time block and also the left boundary of the (i + 1)-th space-time block.

Figure 16 .
Figure 16.Error in the pressure ϵp sampled at x = 0.5 for (a) various values of m with fixed values of ξ 2 = 0.05, and (b) different ξ 2 (m = 100 is fixed in this test).

Figure 17 .
Figure 17.(a) The geometry of the domain and the boundary conditions for the synthetic jet simulation.(b) A snapshot of the vorticity field (non-dimensionalized by multiplying with the length scale of the orifice and dividing by maximum velocity Umax) of the 2D synthetic jet flow.
field and error at the center of the domain (i.e.(x * , y * ) = (0.25, 0.25)) over time.The non-dimensional pressure (p * ) and error (ϵ * p ) calculated frame by frame by a 2D pressure Poisson solver, sampled at the same location, are also shown in figure 20 for comparison.It can be observed from figures 20(a) and (b) that the reconstructed pressure and corresponding error are increasingly smoother with larger ξ 2 .

Figure 20 .
Figure 20.The pressure (a) and the corresponding error in the pressure (b) sampled at (x * , y * ) = (0.25, 0.25).The insets show the zoomed-in views of a portion of the samples.

Figure 23 .
Figure 23.Sample of pressure (a) and the corresponding error in the pressure (b) at (x * , y * ) = (0.25, 0.25).Zoomed-in views of a portion of the samples are shown in the insets.

Figure 25 .
Figure 25.Ar * as function of ξ 2 (a) and T (b) for various combination of values of X, T, F, and ξ 2 with different boundary conditions.The dashed lines represent the cases where mixed boundary conditions are applied, and the solid lines indicate the domain with pure Dirichlet boundary conditions.

Figure 26 .
Figure 26.Ar * as a function of ξ 2 for various combinations of values of X, Y, T and F. The solid lines (k = 1) and dashed lines (k = 3) represent the theoretical upper bound based on inequality (23).The square symbols are based on the numerical tests with k = 3 illustrating the error amplification ratio when the ξ 2 is varied.

Figure 27 .
Figure 27.Ar * as function of ξ 2 for various combinations of X, Y, T, and F. The solid lines represent the (1 + 1)D theoretical upper bound based on inequality (24).The dash lines represent the (2 + 1)D theoretical upper bound based on inequality (35).