Local pressure-correction for the Navier-Stokes equations

This article presents a novel local pressure correction method for incompressible fluid flows and documents a numerical study of this method. Pressure correction methods decouple the velocity and pressure components of the time-dependent Navier-Stokes equations and lead to a sequence of elliptic partial differential equations for both components instead of a saddle point problem. In some situations, the equations for the velocity components are solved explicitly (with time step restrictions) and thus the elliptic pressure problem remains to be the most expensive step. Here, we employ a multiscale procedure for the solution of the Poisson problem related to pressure. The procedure replaces the global Poisson problem by local Poisson problems on subregions. We propose a new Robin-type boundary condition design for the local Poisson problems, which contains a coarse approximation of the global Poisson problem. Accordingly, no further communication between subregions is necessary and the method is perfectly adapted for parallel computations. Numerical experiments regarding a known analytical solution and flow around cylinder benchmarks show the effectivity of this new local pressure correction method.

multigrid is often the best choice for the pressure correction methods.Issues may arise in the context of highly parallel architectures because of the limits of scaling for large numbers of cores.
Lately, several authors 6,7 aimed to solve this Poisson problem in a more computationally efficient way by modifying the finite dimensional approximation spaces.Recently, a local pressure correction method 8 for the Stokes system was proposed.The method employs a domain decomposition technique for the Poisson problem and opens up the possiblity to solve Neumann type Poisson problems in subregions.Here, we numerically investigate the case when Robin boundary conditions are employed for the local Poisson problems.First, we solve the global Poisson problem approximatively on a coarse grid.Consequently, on the same grid with velocity components, solutions to local Poisson problems are approximated .The coarse solution and its gradients enter in the boundary conditions of local problems.This choice of boundary conditions is motivated by the domain decomposition method of Lions, 9 with the adjustment that replaces solutions from neighbour subregions by coarse global solutions in order to prevent iterations.Therefore, parallel computation of the local pressure corrections is very natural and straightforward.Interested reader may consult to Knopp et al. 10 for an extensive review of domain decomposition methods for incompressible flows.
At this point, we would like to emphasize that pressure correction methods are not the most accurate schemes to solve flow problems, in general.Solving the equations for pressure and velocity in a coupled fashion leads to more accuracy, but also to higher numerical costs. 11Furthermore, recently so-called pressure-robust methods 12 has gained popularity, as they lead to a decoupling of the pressure error from the velocity error.This may be beneficial for flows with strong Coriolis forces.We stress the fact that the design of local pressure correction schemes is motivated by durable popularity of pressure correction methods in the literature and in many commercial CFD codes as well.
In Section 2, we introduce the local pressure-correction method with Robin boundary conditions.The finite element discretization and fully discrete solution procedure for the local pressure-correction method is given in Section 3. The numerical studies are presented in Section 4. We finalize the article with a conclusion and outlook.

LOCAL PRESSURE-CORRECTION FOR THE NAVIER-STOKES EQUATIONS
We consider the time-dependent incompressible Navier-Stokes equations in a Lipschitz domain Ω ⊂ R d , d ∈ {2, 3} with homogeneous Dirichlet conditions.In terms of velocity field u ∶ Ω → R d , and pressure p ∶ Ω → R, the equations read where  > 0 is a positive, constant viscosity and f ∶ Ω → R d is a forcing term.The weak solutions to problem (1-4) are sought in the space Many different pressure correction methods can be found in the literature, starting with the pioneering work of Chorin 1 and Temam. 2 Instead of reporting here on the different pressure correction schemes, we refer to the overview article of Guermond et al. 3 In these methods, each time step starts with a computation of a not necessarily divergence-free predictor velocity field, followed by a Poisson problem for the pressure (or a pressure update), and terminates with a projection of the previously computed velocity onto a divergence free one.Although the predictor velocity often demands an implicit method for stability reasons, in certain applications (e.g., in climate research) with the need of (relatively) small time steps due to accuracy reasons, the predictor step can be formulated in an explicit way.In this case, the Poisson problem for the pressure (update) is relatively expensive due to the need of an implicit solver with an ill-conditioned stiffness matrix.The condition number  becomes worse on finer meshes:  ∼ h −2 for mesh sizes of order h.Therefore, parallelization on multi-core machines or parallel computers may help to reduce the simulation time but always require suitable further iteration techniques to account for the elliptic character of the Poisson problem.Recently, a local pressure correction method 8 was proposed.The idea is basically to solve local Poisson problems, independently from each other, with boundary values obtained from a global coarse-grid problem.It was shown that for the Stokes problem the same order of convergence as the original scheme can be obtained.While the boundary conditions for the local problems are of Neumann-type, we propose in this work Robin-type boundary conditions for the subdomain problems.Although the following domain decomposition procedure is applicable to all existing pressure correction methods in the literature, we restrict ourselves to the rotational incremental pressure correction proposed by Timmermans et al. 13 which is the most accurate one beyond others.
We consider a non-overlapping decomposition of Ω: We denote S ∶= H 1 (Ω) and the local pressure corrections in each subdomain Ω i will be sought in the space

Algorithm 1. The local pressure correction
The local pressure correction starts with a given u 0 and p 0 and for m ≥ 1: Step 1: Compute the velocity predictor ũm from momentum equation with a BDF(2)-scheme and a given pressure Step 2a: Find qm ∈ S such that qm ⋅ n = 0 on Ω, ( qm , 1) Ω = 0.
Step 2b: Compute local pressure corrections q m i from Poisson equations on sub-domains where  > 0 is an algorithmic parameter to be defined later.
Step 3: Pressure update: Step 4: Velocity correction: Step 5: If m = M stop.Otherwise, increase m by one, increase t by k, and go to Step 1.
Remark 1.If N = 1 we recover the original scheme. 13emark 2. The Step 2a of Algorithm 1 consists of calculation of the global correction, which enters in the Robin boundary conditions of the local Poisson problems in Step2b.Note that Step 2a does not involve an equality in the pointwise sense, but in an (coarse) approximation sense.Therefore, we use the symbol "≈" in Step 2a.
Remark 3. In case of a pointwise equality in Step 2a and q i = q m in Step 2b, Algorithm 1 corresponds to the original scheme for any N ∈ N which has the following accuracy for sufficiently smooth u and p: While the errors of predictor velocity and corrected velocity in l 2 (0, T; L 2 ) norm at discrete time points are of order k 2 , in l 2 (0, T; H 1 ) norm they are of order k 3/2 .The pressure error in l 2 (0, T; L 2 ) norm is of order k 3/2 as well.For a detailed proof, we refer to Guermond and Shen. 14mark 4. For building the sum in Equation ( 5), the individual q m i , defined on Ω i , are extended to Ω by setting it to zero on Ω ⧵ Ω i .Accordingly, the gradient ∇q m in ( 7) is only defined on ⋃ i Ω i but not on interior boundaries Ω i ∩ Ω j for i ≠ j.Therefore, this velocity update has to interpreted in weak form and by integration by parts in order to shift the gradient onto the test function, see the next section.This method can be considered as a domain decomposition method without iteration within one time-step: The global Poisson problem in Step 2a will be approximated by solving for a discrete pressure update on a coarse mesh.Afterward, the local problems of Step 2b can be solved in parallel without any communication between the subdomains.Similar to other domain decomposition methods (see, e.g., Bjorstad et al. 15 , Chan 16 and references therein), the kind of boundary conditions for the subdomains matters.Here, we propose Robin-type conditions with a right hand side obtained from the global problem in Step 2a.Due to the absence of communication between the sub-domains, we expect a speed up for parallel computations.

FINITE ELEMENT DISCRETIZATION
Let  h (Ω) be a shape-regular, admissible decomposition of Ω into simplices.
Let K be a reference element of  h (Ω) and P k ( K) the space of all polynomials on K with maximal degree k in each coordinate direction.We employ the continuous conforming finite element space Discrete subspace pair for velocity and coarse pressure corrections are Moreover, discrete local pressure corrections are sought in S i,h ∶= P  h (Ω i ),k−1 ∩ S i .After employing the finite dimensional spaces introduced above, in Algorithm 2, we present the spatially discrete version of the Algorithm 1.
Remark 5. Note that, the computation of u m h in the Step 4 of Algorithm 1 was eliminated in Algorithm 2 by substituting it in (8).Moreover, since p m h , q m−1 h , q m−2 h ∈ Q h integration by parts was performed over pressure and pressure updates in Equation (8).The interface term ,  ⋅ n i ) Γ i is neglected, since the mean values of the variables involved in this term are arbitrary on each subdomain.

NUMERICAL EXAMPLES
In order to compare the accuracy of the proposed method numerically, we consider two 2d examples from the literature.First one is with a known analytical solution in terms of trigonometric functions and the second one is the so-called flow around cylinder benchmark problem where lift and drag forces, and as well as certain pressure changes are quantities of interest.We choose  0 = 10 if not stated otherwise.All experiments were performed on FreeFem++. 17The local pressure correction starts with a given u 0 h and p 0 h and for m ≥ 1: Step 2a: Find q m H ∈ S H such that Step 2b: For each i ∈ {1, … , N} find q m i,h ∈ S i,h such that where  =  0 h −1 and  0 ∼ (1).
Step 3: Find q m h and p m h ∈ Q h such that Step 4: If m = M stop.Otherwise, increase m by one, increase t by k, and go to Step 1.
We chose constant time step sizes k = 0.1 ⋅ 2 −n , n ∈ {0, … , 8} and use Taylor-Hood P 2 ∕P 1 finite elements on a triangulation with mesh size h = 2 −7 .The elliptic coarse grid problem for q H is carried out by P 1 -elements on a triangulation with mesh size H ∈ {8h, 16h}.Therefore, solving this problem is a factor 64 or 256 smaller (with respect to the number of unknowns) than the one for the pressure update of the original pressure correction scheme.That is to say, the computation of q H is basically for free.For this example, we will study the velocity error in the following norm and semi-norm The pressure error is considered in the following semi-norm We note that there is no way to really determine the mean value of pressure on each subdomain, which would correspond to the one of the global pressure.In the following two subsections, we present the results for the diffusion dominant case ( = 1) and for the convection dominant case ( = 10 −3 ), respectively.In Section 4. on the accuracy of the local pressure correction scheme.Afterward, in Section 4.1.5,results with the Neumann type local pressure corrections will be given.

Results for the diffusion dominant case
In the first test, we consider the Navier-Stokes problem with viscosity  = 1 for a variety of number of subdomains N ∈ {1, 4, 16, 64, 256}, where N = 1 corresponds to the original scheme without subdomains.In Figures 1 and 2, we present convergence rates of ||u in dependence of the time step size k.It turns out that there is no difference in the accuracy of local pressure correction method when compared to the original scheme.Moreover, we note that the convergence rate in norm ||u − ũh || l 2 (0,T;H 1 (Ω)) is (k 2 ), while the theory predicts (k 3∕2 ).The convergence rates in norms ||u − u h || l 2 (0,T;L 2 (Ω)) and ||u − ũh || l 2 (0,T;L 2 (Ω)) are (k 2 ) as predicted by the theory for original scheme.Furthermore, the pressure error with respect to the spatial gradient ||p − p h || l 2 (0,T;H 1 (Ω)) for the local pressure correction leads to the same accuracy as the original scheme, independent of the number of subdomains.We observe that the curve flattens for smallest time steps as a consequence of the dominating spatial error.

Results for the convection dominant case
As next, we modify the viscosity to  = 10 −3 .It is well known that for smaller values of viscosity, the problem (1) becomes convection dominant and numerical approximation becomes more challenging.Moreover, a different scaling of the diffusive term leads to a different scaling between pressure and velocity.The error results are presented in Figures 3 and 4.
For H = 8h, again, there is no difference in the accuracy of the local pressure correction method and the original scheme.The convergence rates are the same as in the case  = 1.On the other hand, for H = 16h, we observe that, only for small time steps, there is an accuracy loss in the norms ||u − u h || l 2 (0,T;L 2 (Ω)) and ||u − ũh || l 2 (0,T;H 1 (Ω)) .This behavior suggests that the coarse grid problem cannot be solved on triangulations with an arbitrarily large mesh size H.

Results for the case of dominating spatial discretization error
Now we choose  = 10 −3 and h = 2 −5 .In Figure 5, we present a comparison of the original pressure correction scheme (N = 1) and the local pressure correction scheme with N = 16 and H = 8h.We observe that the spatial pressure gradients from both schemes stagnate at k = 0.025.Both methods lead to same errors for predictor velocity in the norm ||u − ũh || l 2 (0,T;L 2 (Ω)) .For the smallest time step size the errors of the local pressure correction method (N = 16) in the norms ||u − u h || l 2 (0,T;L 2 (Ω)) and ||u − ũh || l 2 (0,T;H 1 (Ω)) are the same of the original pressure correction scheme.

The choice of parameter 𝝉 0
Now we consider different choices for the parameter  0 .We choose N = 64, h = 2 −7 and H = 16h.From Section 4.1.2,it is clear that this is a challenging setting for the local pressure corrections.In Figure 6, we present the errors for  0 ∈ {10, 10 4 , 10 8 } and compare to the original pressure correction method (N = 1).It becomes evident that the choice of  0 is of small significance for k ≥ 0.0125.However, as the time step size chosen smaller, the error of corrector velocities in the norms ||u − ũh || l 2 (0,T;L 2 (Ω)) and ||u − ũh || l 2 (0,T;H 1 (Ω)) grows faster for  0 ∈ {10 4 , 10 8 } when compared to  0 = 10 or the original pressure correction method (N = 1).On the other hand, from our experiences  0 < 10 leads to an unstable scheme, in general.A theoretical analysis for the design of  0 remains open.

Comparison with local pressure corrections of Neumann type
As next we present results obtained with local pressure corrections of Neumann type. 8In this formulation, the Step 2b reads: for each i ∈ {1, … , N} find q m i,h ∈ S i,h such that where The results depicted in Figure 7 for h = 2 −7 , H = 8h, and N = 64 show that the local pressure corrections of Neumann type yield to a large loss of accuracy for smaller time step sizes.Moreover, a comparison of Figure 3(N = 64) and Figure 7indicates that, for k ≥ 0.025, the errors obtained from both local pressure corrections of Neumann type and Robin type are in agreement.However, for smaller time step sizes, the local pressure corrections of Robin type outperforms the Neumann type as they exhibit the same asymptotic behavior with the original scheme.Therefore, we deduce that the Robin type local correction projections are favorable when compared to the Neumann type.

Benchmark problem "Flow around cylinder"
In this subsection, we present our numerical results concerning the so-called flow around cylinder benchmark problem according to the domain configuration given in Schäfer et al. 11 and compare with the reference results in John 18 .For interested readers, we refer to John et al. 19 for an overview of time integration techniques related to this problem.Domain of interest is Ω = [0, 2.2] × [0, 0.41] ⧵ S where S is the circle around the point (0.2, 0.2) with radius 0.05.Moreover, we set I : = (0, T] with T = 8 and f = 0.The inlet and outlet boundaries are denoted by Γ in = {(0, y) ∈ R 2 ; 0 ≤ y ≤ 0.41} and Γ out = {(2.2,y) ∈ R 2 ; 0 ≤ y ≤ 0.41}, respectively.Boundary conditions are given as The mean inflow velocity is U(t) = sin(t∕8) and thus U max = 1.Therefore, based on L = 0.1 (the diameter of S) and  = 10 −3 , the Reynolds number of the flow is found to be 0 ≤ Re ≤ 100.Drag and lift cofficients at the cylinder are calculated with the formula  where  = 1 (fluid density) is a known parameter.Here u t S denotes the tangential velocity, n = (n x , n y ) T is the normal vector on S directing into Ω and t S = (n y , − n x ) T is the tangential vector.Moreover, the pressure difference between the front and the back of the cylinder D p (t) = p(t; 0.15, 0.2) − p(t; 0.25, 0.2), is another parameter of interest.The domain Ω is decomposed into 24 subdomains as illustrated in Figure 8. Two spatial discretization settings are considered whereas both employ the same coarse mesh  H (Ω) (see Figure 9).While H = 4h for Mesh 1, H = 8h was chosen for Mesh 2. Computation of the coarse grid correction q H is then negligible regarding the numerical costs, when compared with the original correction.The local Poisson problems on each subdomain are solved on  h (Ω i ).The degrees of freedom corresponding to each local pressure correction q m i,h ∈ S i,h are illustrated in Figure 10.Note that the choice of subdomains is rather flexible and certain modifications can decrease the size of local problems.In Table 1, we present the number of degrees of freedom for velocity ũm h ∈ V h , original (global) pressure correction q m h ∈ Q h , and the coarse pressure correction q m H ∈ S H .A challenging aspect of this setting is: there are four subdomains around the cylinder, which means the pressure differences are determined on four independent local problems.In Figures 11 and 12, we present drag and lift coefficients and pressure differences computed on Mesh 2. For local pressure corrections, average pressure differences on the subdomains (5, 6) and (7, 8) are considered.Note that moderate differences are only present in the lift coefficient.In Figure 11, we zoom around the region where the lift coefficient attains its maximum value.The reference maximum value of the lift coefficient 18 is depicted by a point.We observe that, for large time-steps, there is a moderate difference between lift coefficients of local pressure and original pressure correction methods.However, this difference becomes less when the time-step is chosen smaller and both methods converge to the reference value.This behavior can be observed in Table 2 in more detail for both mesh settings.Regarding the drag coefficient, we stress the fact that the original pressure correction method does not converge to the reference value, but leads to a relative error of approximately 2 ⋅ 10 −4 on Mesh 2. Nevertheless, the drag coefficient and pressure differences resulting from local pressure correction method converges to the one from original pressure correction method on Mesh 1 and 2.

CONCLUSIONS
We presented a numerical study of the local pressure correction method for approximating the solutions to Navier-Stokes equations.The method is based on local projections of the velocity onto divergence-free velocities.The global Poisson problems are replaced with a (very) coarse global and several local Poisson problems.The resulting local Poisson problem can be solved completely in parallel.We numerically investigated that for velocity the discretization error of the method has the same asymptotic behavior as the original pressure correction scheme.Moreover, we presented our computational results regarding flow around cylinder benchmark, where physical quantities such as lift and drag forces and pressure differences are of interest.The results show that the local pressure correction scheme is qualitatively very promising.One of the next steps will be to solve the subproblems in parallel and to document on the speed up of the new scheme.Moreover, a theoretical error analysis of the scheme is inevitable as it would lead to a robust design of the parameter  0 .

Algorithm 2 .
The local pressure correction, fully discrete case

F I G U R E 1
1.3, we present results for the case when spatial discretization error dominates.Later on in Section 4.1.4,we present the results for the effect of different  0 Numerical results for  = 1, H = 8h (Section 4.1.1)[Color figure can be viewed at wileyonlinelibrary.com]

F I G U R E 2
Numerical results for  = 1, H = 16h (Section 4.1.1)[Color figure can be viewed at wileyonlinelibrary.com]

F I G U R E 3
Numerical results for  = 10 −3 , H = 8h (Section 4.1.2)[Color figure can be viewed at wileyonlinelibrary.com]

F I G U R E 4
Numerical results for  = 10 −3 , H = 16h (Section 4.1.2)[Color figure can be viewed at wileyonlinelibrary.com]

F I G U R E 7 F
figure can be viewed at wileyonlinelibrary.com]

F I G U R E 11
Mesh 2: Obtained lift coefficients (left) for the original (Global) scheme and the novel (Local) scheme with different time steps k = 0.01 ⋅ 2 −s .Zoom around maximal lift coefficients (right) [Color figure can be viewed at wileyonlinelibrary.com]F I G U R E 12 Mesh 2: Obtained drag coefficients (left) and pressure differences (right) with different time steps k = 0.01 ⋅ 2 −s from the original and the novel scheme.All curves are nearly not distinguishable [Color figure can be viewed at wileyonlinelibrary.com]TA B L E 2 Maximal lift and drag coefficientes and pressure differences; Global scheme (N = 1, in bold) and Local scheme (N = 24, plain) … , N}.The interfaces are denoted by Γ i ∶= Ω i ⧵ Ω and the outer normal on the boundary Ω i by n i .For the sake of simplicity, we take constant time step size k = t m − t m − 1 > 0 for all m, where 0 = t 0 , t 1 , … , t M are M + 1 discrete time points.Hence t m : = mk for 0 ≤ m ≤ M = T/k.The time discrete solution at time t m will be denoted by (u m , p m ).In general, pressure correction schemes solve in each step a Poisson problem on Ω equipped with homogeneous Neumann boundary conditions.Consideration of several non-overlapping subdomains Ω i , i ∈ {1, … , N} instead of Ω brings out the issue of determining consistent boundary conditions for local Poisson problems on subdomains.Our approach is to equip local Poisson problems with Robin type boundary conditions, where a computationally cost efficient approximation of the Poisson problem on Ω and its gradients enter on the right hand side of the boundary data.The details of the resulting local pressure correction method are presented in Algorithm 1.
Number of degrees of freedom