Parametric solutions of turbulent incompressible flows in OpenFOAM via the proper generalised decomposition

An a priori reduced order method based on the proper generalised decomposition (PGD) is proposed to compute parametric solutions involving turbulent incompressible flows of interest in an industrial context, using OpenFOAM. The PGD framework is applied for the first time to the incompressible Navier-Stokes equations in the turbulent regime, to compute a generalised solution for velocity, pressure and turbulent viscosity, explicitly depending on the design parameters of the problem. In order to simulate flows of industrial interest, a minimally intrusive implementation based on OpenFOAM SIMPLE algorithm applied to the Reynolds-averaged Navier-Stokes equations with the Spalart-Allmaras turbulence model is devised. The resulting PGD strategy is applied to parametric flow control problems and achieves both qualitative and quantitative agreement with the full order OpenFOAM solution for convection-dominated fully-developed turbulent incompressible flows, with Reynolds number up to one million.


Introduction
Parametric studies involving flows of industrial interest require robust computational fluid dynamics (CFD) solvers and efficient strategies to simulate multiple queries of the same problem.
When the simulation requires testing a large number of different configurations -e.g. for shape optimisation, uncertainty quantification, inverse and control problems -numerical strategies to reduce the cost of the overall computation are critical. Reduced order models (ROM) [15,16] construct an approximation of the solution in a lower dimensional space, for which an appropriate basis needs to be devised.
It is known that numerical difficulties arise when reduced basis (RB) and proper orthogonal decomposition (POD) techniques are applied to convection-dominated problems [17][18][19]. This is especially critical in the context of flow simulations when the Reynolds number is increased and turbulent phenomena need to be accounted for. More precisely, the most relevant POD modes are associated with the highest energy scales of the problem under analysis, whereas small scales, which play a critical role in the dissipation of turbulent kinetic energy, are poorly represented by POD-ROM [20]. To remedy this issue, closure models stemming from traditional description of turbulence have been extended to ROMs, leading to Galerkin projection-based POD with dynamic subgrid-scale [21], variational multiscale [22,23], k − ω SST [24] models and to a certified Smagorinsky RB strategy [25]. Moreover, strategies to improve efficiency and accuracy of POD-ROM in the context of realistic and turbulent flows have been proposed by coupling the projectionbased framework with residual minimisation [26], nonlinear least-squares optimisation [27], interpolation based on radial basis functions [28] and a constrained greedy approach [29]. In the context of machine learning-based reduced order models [30][31][32], a strategy coupling a traditional projection-based POD for velocity and pressure with a data-driven technique for the eddy viscosity has been recently proposed in [33].
All above contributions involve the development of a posteriori ROMs, namely RB and POD, in which the basis of the low-dimensional approximation space is computed starting from a set of snapshots. On the contrary, PGD [34,35] constructs a reduced basis of separable functions explicitly depending on space and on user-defined parameters, with no a priori knowledge of the solution of the problem. The resulting PGD computational vademecum provides a generalised solution which can thus be efficiently evaluated in the online phase via interpolation in the parametric space, that is, no extra problem needs to be solved in the low-dimensional reduced space as in POD. In the context of flow problems, PGD was originally utilised to develop efficient solvers for the incompressible Navier-Stokes equations by separating spatial directions [36] and space and time [37,38]. In addition, problems involving parametrised geometries have been solved using PGD [39,40], with special emphasis on incompressible flows in geometrically parametrised domains [41][42][43][44]. To foster the application of a priori model order reduction techniques to problems of industrial interest, a non-intrusive PGD implementation in the CFD software OpenFOAM has been recently proposed in [45] to solve parametrised incompressible Navier-Stokes equations in the laminar regime.
Following the work on PGD for convection phenomena [46,47] and for viscous incompressible Navier-Stokes flows [45], the present contribution proposes the first a priori ROM for turbulent incompressible flows. Similarly to the previously cited contributions on a posteriori ROMs [21][22][23][24][25][26][27][28][29][30][31][32][33], a priori ROMs can also be devised using a variety of techniques to model turbulence. The proposed PGD strategy considers the one-equation SA turbulence model and it constructs a separated representation of velocity, pressure and eddy viscosity to solve the RANS-SA equations in OpenFOAM. More precisely, the PGD-ROM methodology mimics the structure of the simpleFoam algorithm with SA turbulence model, resulting in a minimally intrusive approach within OpenFOAM. The resulting strategy thus provides a generalised expression of the velocity, pressure and eddy viscosity fields, explicitly depending on user-defined parameters, for convection-dominated incompressible flows. Alternative PGD-ROM strategies may thus be obtained by substituting the proposed separated form of the SA equation with an appropriate PGD solver for the selected turbulence model.
The remainder of this paper is organised as follows. Section 2 recalls the full order RANS-SA equations and the corresponding cell-centred FV approximation utilised by OpenFOAM. The rationale of the PGD-ROM for the turbulent incompressible Navier-Stokes equations is introduced in section 3, where the details of the algorithms to devise the separated velocity-pressure approximation of the flow equations (PGD-NS) and the separated form of the eddy (PGD-SA) and turbulent (PGD-ν t ) viscosities via the SA equation are presented. Numerical experiments involving flow control in external aerodynamics, in two and three dimensions, with Reynolds number ranging up to 1, 000, 000 are reported in section 4. Finally, section 5 summarises the contributions of this work and two appendices report additional technical details on the employed PGD algorithm and on the expressions of the coefficients appearing in the spatial and parametric iterations of the alternating direction scheme for the PGD solvers of the RANS and the SA equations.

The Reynolds-averaged Navier-Stokes equations and the Spalart-Allmaras turbulence model
To simulate turbulent incompressible flows using the RANS equations, the velocity-pressure pair (u, p) is decomposed into a mean flow component (U , P ) and a perturbation (u , p ), that is u=U +u and p=P +p . Given an open bounded computational domain Ω ⊂ R d in d spatial dimensions, the boundary ∂Ω is partitioned such that ∂Ω=Γ in ∪ Γ w ∪ Γ out , where the three disjoint portions Γ in , Γ w and Γ out denote inlet surfaces, material walls and outlet surfaces, respectively. The steady-state RANS equations for the mean flow variables (U , P ) are given by where I d denotes the d×d identity matrix, ν represents the physical viscosity of the fluid and ν t is the turbulent viscosity introduced in the momentum equation to model the perturbations to the mean flow due to turbulence. The boundary conditions for the flow equations (1) impose the velocity profile U in on the inlet surface Γ in , no-slip Dirichlet data on fixed material walls Γ w and homogeneous Neumann data on the outlet surface Γ out .
In order to describe the turbulent viscosity ν t , the one-equation SA turbulence model [13] introduces the relation where ν is the eddy viscosity and f v1 is an appropriately defined spatial function, see e.g. [13,48,49], reported in equation (4). Under the assumption of fully turbulent flows, the trip term controlling the transition between laminar and turbulent regimes in the SA model is neglected and the eddy viscosity ν is obtained as the solution of where ν D is the eddy viscosity Dirichlet datum, d represents the distance of a given point in the domain from the closest physical wall and c b1 , c b2 , c w1 and σ are four scalar constants.
Moreover, S and f w are the spatial functions associated with the production and the destruction of eddy viscosity, respectively, whereas the operators on the left-hand side of equation (3) model convection, diffusion and cross-diffusion phenomena [13].

A finite volume formulation of the RANS-SA equations
In order to discretise the turbulent Navier-Stokes equations, OpenFOAM cell-centred finite volume rationale is considered [14]. The computational domain is subdivided in N cells where (U , P ) are cell-by-cell constant approximations of the velocity and pressure fields, respectively, and U =U in on Γ in and U =0 on Γ w .
In a similar fashion, the cell-centred finite volume approximation of the SA equation (3) is: compute ν constant in each cell such that ν = ν D on Γ in ∪ Γ w and it holds

A turbulent Navier-Stokes solver in OpenFOAM
OpenFOAM strategy to solve the RANS equation with SA turbulence model relies on a staggered approach. First, the flow equations (5) are solved using a seed value of ν t . More precisely, the integrals over each cell in (5) are approximated by means of the corresponding fluxes across the boundaries of the cell [4,5]. In addition, the semi-implicit method for pressure linked equations (SIMPLE) algorithm [50], that is, a fractional-step Chorin-Temam projection method [51,Sect. 6.7], is utilised to handle incompressibility. Also, a relaxation approach is employed for the nonlinear convection term. Second, the velocity field U obtained using simpleFoam is employed to compute the quantities in (4) and to solve the SA equation (6). It is worth noting that equation (6) is highly nonlinear and a relaxation strategy is also required in OpenFOAM to improve the convergence of the numerical algorithm [14]. Finally, the updated value of the turbulent viscosity ν t is determined according to equation (2) and the simpleFoam routine is utilised to recompute the turbulent velocity and pressure fields.

Proper generalised decomposition for parametric turbulent flow problems
In the context of parametric studies, viscosity coefficient, reference velocity or boundary conditions of the problems may depend on a set of M user-defined parameters µ=(µ 1 , . . . , µ M ) T . The solution of the RANS-SA equations is thus denoted by the velocity-pressure pair (U (x, µ), P (x, µ)) and the eddy viscosity ν(x, µ), which are now functions of the spatial, x ∈ Ω ⊂ R d , and parametric, µ ∈ I ⊂ R M , variables. More precisely, U (x, µ), P (x, µ) and ν(x, µ) fulfil the high-dimensional RANS-SA equations obtained by integrating (5)- (6) in the parametric space I.
The PGD-ROM strategy described in this section relies on the construction of a separated approximation (U n PGD , P n PGD ) of the velocity and pressure fields and a separated representation of any additional variable introduced by the employed turbulence model. The computation of the former is performed via a PGD solver for the incompressible Navier-Stokes equations, see section 3.3, whereas for the latter a separated formulation of the involved turbulence equations is required.
Considering the one-equation SA turbulence model introduced in section 2, a separated representation ν m PGD of the eddy viscosity is obtained from a PGD solver of the SA equation as described in section 3.4. In addition, a separated representation ν q t,PGD of the turbulent viscosity is devised in section 3.5 exploiting its relation with the eddy viscosity. It is worth noticing that this framework is general and could be adapted to other descriptions of turbulence by devising separated representations of the involved variables via appropriately defined PGD solvers of the equations in the turbulence model.
The global set of variables involved in the PGD approximation of the RANS-SA equations is thus (U n PGD , P n PGD ), ν m PGD and ν q t,PGD . As classical in PGD [34], each variable is constructed as a sum of separable modes, each being the product of functions that depend on either the spatial or one of the parametric variables µ j , j = 1, . . . , M . For the sake of simplicity, only space, x, and parameters, µ, are henceforth separated. It is worth noticing that the final number of modes needed for the PGD approximations, denoted by the super-indexes n for the velocity and the pressure, m for the eddy viscosity and q for the turbulent viscosity, is not known a priori and, in general, it is different for each of the involved variables. More precisely, the number of terms in the PGD expansion is automatically determined by the algorithm which stops the enrichment procedure when a user-defined stopping criterion is fulfilled [34]. Classical definitions of this stopping criterion include the relative amplitude of the last computed mode with respect to the first one or to the sum of all previously computed terms [45].

Separated representation of the flow and the turbulent variables
First, the rank-n separated representation (U n PGD , P n PGD ) of the flow variables is introduced. Following [45], the computation of each PGD mode is split into a prediction and a correction step. More precisely, the PGD approximation for the flow variables is defined as where U n−1 PGD and P n−1 PGD feature the contributions of the previous n−1 PGD modes, σ n U f n U φ n and σ n P f n P φ n represent the predictions of the n-th mode and σ n U ∆(f n U φ n ) and σ n P ∆(f n P φ n ) are the corresponding correction terms. The coefficients σ n U and σ n P denote the amplitudes of the n-th velocity and pressure mode, respectively. Similarly, the rank-m PGD separated form of the eddy viscosity is given by with ν m−1 PGD containing the previous m−1 terms in the PGD approximation, σ m ν f m ν ψ m and σ m ν ∆(f m ν ψ m ) being the prediction and the correction of the m-th mode, respectively, and σ m ν its amplitude. Both the PGD approximations of the flow variables (7a) and of the eddy viscosity (7b) are devised solving appropriate separated forms of the corresponding equations, as presented in the following sections. More precisely, U n PGD and P n PGD are obtained from a separated form of the Navier-Stokes equations (5), see section 3.3, whereas ν m PGD is determined via a PGD projection of the SA equation (6), as reported in section 3.4.
Finally, a separated representation of the turbulent viscosity is obtained from the relation (2), leading to the rank-q PGD expansion where ν q−1 t,PGD represents the PGD approximation obtained from the previous q−1 terms, f q t and ξ q denote the q-th normalised spatial and parametric modes, respectively, and σ q t its corresponding amplitude. It is worth noticing that, contrary to the flow variables and the eddy viscosity, the computation of the turbulent viscosity ν q t,PGD does not involve the solution of a differential equation and only requires elementary arithmetic operations, whence the predictor-corrector approach is substituted by the classical PGD separation in equation (7c).
Remark 1. Following [41], the same scalar parametric function φ(µ) is selected for both velocity and pressure. On the contrary, different scalar functions ψ(µ) and ξ(µ) are considered for the separated approximations of the eddy and turbulent viscosity, respectively.
The corrections of the PGD modes introduced in equation (7a) and (7b) for the computation of the current mode feature variations ∆ in the spatial and in the parametric functions, namely Of course, higher order contributions (e.g., ∆f U (x)∆φ(µ) for the velocity) could also be considered in the definition of the corrections (8) of the PGD modes. Nonetheless, the importance of these extra terms is negligigle with respect to both the predictions f n U , f n P , φ n , f m ν and ψ m (µ) of the modes and the first-order corrections introduced in equation (8) and they are thus neglected henceforth. Finally, for the sake of readability, the following compact expressions are introduced to define the corrections of the velocity, pressure and eddy viscosity to be computed by the PGD algorithm where the dependence of the spatial and parametric modes on x and µ has been omitted. A comparison of the classical PGD algorithm, see e.g. [34], and the predictor-corrector strategy employed in this work is reported in A.

Remark 2.
Upon convergence of the alternating direction algorithm, the above mentioned predictions and corrections are combined into the n-th spatial and parametric modes denoted by f n U and φ n , with f n U = φ n =1. In this context, the coefficient σ n U encapsulates the information on the amplitude of the mode. From a practical viewpoint, the solution of the alternating direction method is the pair (f U n , φ n ),f U n being the spatial mode before the normalisation procedure. The corresponding normalised spatial mode is thus given by f n U =f U n /σ n U , with σ n U := f U n . The details of the normalisation procedure for the classical and the predictor-corrector PGD strategies are presented in A. For the simulations in section 4, the normalisation procedure has been performed using the L 2 norm.

A minimally intrusive PGD implementation of a parametric solver for turbulent incompressible Navier-Stokes flows in OpenFOAM
The proposed minimally intrusive parametric solver for turbulent incompressible Navier-Stokes flows in OpenFOAM constructs the separated expressions U n PGD , P n PGD , ν m PGD and ν q t,PGD exploiting the numerical discretisation techniques natively implemented in OpenFOAM and validated by the CFD community. For this purpose, the PGD algorithm mimics the segregated structure of simpleFoam with SA turbulence model which involves the following steps [14]: (A) Compute velocity and pressure using a fractional step projection approach to solve (5), given the current approximation of the turbulent viscosity (RANS solver via simpleFoam).
(B) Use the value of the computed velocity to solve (6) and determine the eddy viscosity (SA solver).
(C) Update the turbulent viscosity according to (2) and go to step (A) to restart the computation.
Following this rationale, the corresponding parametric solver is described in algorithm 1. To this end, three routines are introduced to compute U n PGD , P n PGD , ν m PGD and ν q t,PGD : • PGD-NS solves a separated form of the incompressible Navier-Stokes equations to determine the velocity-pressure pair (U n PGD , P n PGD ), given the current PGD expression ν q t,PGD of the turbulent viscosity (Section 3.3); • PGD-SA solves a separated form of the SA turbulence model to determine the eddy viscosity ν m PGD (Section 3.4); • PGD-ν t computes the updated approximation ν q t,PGD of the turbulent viscosity by means of a separated expression of its relation with the eddy viscosity (Section 3.5) and go to PGD-NS to restart the computation.
The minimally intrusive PGD solver in OpenFOAM is thus obtained by integrating the PGD-NS, PGD-SA and PGD-ν t routines in the following computational framework. More precisely, in step (A), the PGD algorithm solves the parametrised flow equations (Algorithm 1 -Line 5) via the non-intrusive PGD-NS strategy, given the current PGD approximation of the turbulent viscosity. The PGD enrichment procedure for the RANS equations continues by alternately solving the spatial and parametric problems until a user-prescribed threshold is achieved by the amplitude of the computed velocity and pressure modes. Once a sufficiently accurate PGD approximation of velocity and pressure is obtained, step (B) computes a separated representation of the eddy viscosity by means of the minimally intrusive PGD-SA block (Algorithm 1 -Line 9). In step (C), the separated representation Algorithm 1 An OpenFOAM implementation of turbulent pgdFoam Require: Stopping criteria η ( = U, P ) and η ν for the PGD enrichment of the flow equations and the turbulence model. Initial accuracy level γ for the PGD enrichment of the turbulent viscosity. 1: Compute boundary condition modes: the spatial mode is solution of (5) using simpleFoam with SA turbulence model and the parametric modes are set equal to 1. 2: Set n ← 1, i ← 0 and initialise the amplitudes σ 1 ← 1. 3: while σ n > η σ 1 do

4:
Set the enrichment threshold η i t for the turbulent viscosity according to equation (10).

5:
Call PGD-NS to compute the spatial, (f n U , f n P ), and parametric, φ n , modes of velocity and pressure and the amplitudes σ n U and σ n P .

11:
end while 12: Call PGD-ν t to compute the spatial, f q t , and parametric, ξ q , modes of the turbulent viscosity and the amplitude σ q t .

14:
Reinitialise the mode counter n ← 1. Reset U n PGD and P n PGD .

17:
end if 18: end while of the turbulent viscosity is recomputed by means of the PGD-ν t routine (Algorithm 1 -Line 12). Finally, the PGD approximation of the flow equations is reset (Algorithm 1 -Line 14) and step (A) restarts the PGD computation of the flow variables using the newly computed separated expression of the turbulent viscosity.
Remark 3. The overall cost of the PGD solver for the parametric turbulent Navier-Stokes equations depends on the number of updates performed to correct the PGD expression of the turbulent viscosity (Algorithm 1 -Lines from 6 to 15). Indeed, the computation of the velocity and pressure modes is restarted from scratch after each update of the turbulent viscosity (Algorithm 1 -Line 14). In this context, computing a large number of modes for velocity and pressure using early approximations of the turbulent viscosity might significantly increase the computational effort of the algorithm with limited accuracy gain. To remedy this issue, effective numerical strategies are devised by introducing an appropriately defined criterion that limits the number of velocity and pressure modes determined in the early stages of the algorithm while allowing a larger number of modes to be computed when the precision of the approximation of the turbulent viscosity increases. More precisely, the PGD enrichment of velocity and pressure for a given expression of the turbulent viscosity stops when the amplitude of the computed modes drops below a userdefined tolerance η i t (Algorithm 1 -Line 6). In addition, this threshold is gradually reduced after each update of the turbulent viscosity (Algorithm 1 -Line 4) allowing to increase the accuracy of the velocity and pressure modes computed by means of the PGD-NS routine, when the PGD-ν t algorithm provides an improved PGD representation of the turbulent viscosity. In the simulations in section 4, the threshold to control the number of velocity and pressure modes computed by the PGD-NS algorithm at the i-th iteration is given by that is, starting from an initial accuracy of 10 −γ , an exponentially decreasing tolerance is defined after each update of the turbulent viscosity. An alternative approach to control the accuracy of the separated representation of eddy and turbulent viscosities may be devised modifying line 6 of algorithm 1 and fixing a priori the number of modes in the PGD approximation of the velocity field required to run PGD-SA and PGD-ν t routines.
In the following sections, the structure of the PGD-NS, PGD-SA and PGD-ν t routines for the computation of (U n PGD , P n PGD ), ν m PGD and ν q t,PGD , respectively, will be detailed.

Proper generalised decomposition of the flow equations
In this section, the spatial and parametric steps of the non-intrusive PGD algorithm applied to the turbulent Navier-Stokes equations (5) are presented. The integral form of the steady parametrised Navier-Stokes equations in the high-dimensional space Ω × I is given by for each cell V i , i=1, . . . , N of the spatial computational domain.
The PGD solver for the incompressible Navier-Stokes equations is obtained by replacing the separated approximations (7a) into equation (11). The first mode (U 0 PGD , P 0 PGD ) is selected to verify the boundary condition on the inlet surface and, more generally, all inhomogeneous Dirichlet boundary data. It is worth recalling that, as classic in PGD approximations [34], the first mode only prescribes the Dirichlet boundary conditions and it does not fulfill the equation under analysis. Then, the following terms of the PGD expansion of velocity and pressure are computed following the greedy strategy described in [45], henceforth named PGD-NS. To compute the n-th modes σ n U f n U φ n + δU n PGD and σ n P f n P φ n + δP n PGD , the highdimensional problem is thus alternatively restricted to the tangent manifold associated with either the spatial or the parametric coordinate. More precisely, a first iteration of the alternating direction scheme is performed to determine the predictions (σ n U f n U φ n , σ n P f n P φ n ) using the classical PGD algorithm [34,52], whereas the following iterations are devoted to compute the separated expressions of the corrections (σ n U δU n PGD , σ n P δP n PGD ) [45].
Remark 4. The segregated algorithm employed by OpenFOAM to simulate turbulent flows requires an initial guess of the turbulent viscosity ν t to solve the flow equations. The turbulent viscosity is later updated after the SA equation is solved. In a similar fashion, the proposed PGD solver utilises the last computed separated approximation ν q t,PGD of the turbulent viscosity to solve the flow equations, before updating it with the result of the PGD-SA and PGD-ν t steps. Thus, the turbulent viscosity ν t in equation (11) is henceforth replaced by its separated approximation ν q t,PGD .
+σ n P f n P φ n denote the estimated rank-n PGD approximations for velocity and pressure obtained as the sum of the converged n−1 terms and the computed predictions of the n-th modes. The unknown increments (σ n U δU n PGD , σ n P δP n PGD ) are determined from equation (11) solving where the unknowns (σ n U δU n PGD , σ n P δP n PGD ) have been gathered on the left-hand side, whereas the right-hand side features the residuals Remark 5. It is straightforward to observe that the equations (12) mimic the structure of the full order Navier-Stokes equations (11). In order to devise a PGD algorithm nonintrusive with respect to the OpenFOAM spatial solver simpleFOAM, in [45] the authors proposed to relax the second and third terms in (12) during the PGD spatial step. More precisely, these terms were evaluated using the last known increment computed during the previous SIMPLE iteration. Hence, they were treated explicitly, as additional terms on the right-hand side of the momentum equation. This approach is also followed in this work, as detailed in section 3.3.1.
An efficient implementation of this strategy can be devised by exploiting the affine decomposition of the forms in (12) and (13), see [53,54], and the separated structure of the unknowns (7a). Nonetheless, it is worth noticing that when the number of modes in the separated representation grows, the overall cost of the PGD strategy may significantly increase. In this context, several variations of the above approach have been proposed in the literature, e.g. by performing a PGD compression step [55,56] to eliminate the redundant information in the previously computed modes or by introducing an update step [57,58] or an Arnoldi-type iteration [59,60] to exploit the intermediate solutions in the alternating direction algorithm.
In the following sections, the restriction of the high-dimensional problem (12) to the tangent manifold associated with the spatial and the parametric coordinates is described. More precisely, at each iteration of the alternating direction algorithm, the spatial step is solved using simpleFoam, whereas the parametric step leads to an algebraic problem solved using a collocation approach.

PGD-NS: the spatial iteration
In order to construct a PGD approximation of the flow equations (12), a separated representation of the data is required [40]. More precisely, the form ν(x, µ)=D(x)ζ(µ) is assumed for the physical viscosity, whereas the PGD approximation The high-dimensional problem (12) is thus restricted to the spatial direction multiplying it by φ n . In addition, setting the value of the parametric function φ n , the increments (σ n U δU n PGD , σ n P δP n PGD ) reduce to φ n (σ n U ∆f U , σ n P ∆f P ), see equation (9). The spatial increments (σ n U ∆f U , σ n P ∆f P ) are thus computed as the FV solution of the flow equations where R n U and R n P denote the spatial residuals of the momentum and continuity equations, respectively, and the coefficients α k , k=0, . . . , 3 and α 7 , reported in B, only depend on user-defined data and parametric functions and can thus be efficiently precomputed.
Following remark 5, an implementation of the PGD spatial solver for the flow equations non-intrusive with respect to the OpenFOAM SIMPLE algorithm is obtained by relaxing the two linear contributions arising from the nonlinear convection term [45]. More precisely, the last two integrals on the right-hand side of the momentum equation in (14) are evaluated using the last increment σ k−1 U ∆f k−1 U computed in the SIMPLE iterations. It is straightforward to observe that the resulting structure of the left-hand side of equations (14) mimics the traditional Navier-Stokes equations (5), whence the PGD spatial iteration is solved using the simpleFoam algorithm, natively implemented in OpenFOAM [45].
The spatial residuals R n U and R n P on the right-hand side of equations (14) are determined starting from the previous PGD terms (U n−1 PGD , P n−1 PGD ) and from the predictions (σ n U f n U φ n , σ n P f n P φ n ) of the n-th mode currently computed, namely It is worth recalling that the factor φ n in the expressions of R n U and R n P follows from the restriction of the residuals (13) defined in the high-dimensional space Ω×I to the tangent manifold associated with the spatial direction [45]. In order to perform an efficient computation of such residuals, the separated expressions of (U n−1 PGD , P n−1 PGD ) as a product of spatial and parametric functions are exploited, leading to where the coefficients α k , k=4, . . . , 6 and α 8 encapsulate the information of the previously computed parametric modes and are defined in B.

PGD-NS: the parametric iteration
In the parametric iteration of the PGD solver for the flow equations, the value of the spatial functions (σ n U f n U , σ n P f n P )←(σ n U [f n U +∆f U ], σ n P [f n P +∆f P ]) is updated and fixed. It follows that the increments associated with the momentum and continuity equations in the parametric step are σ n U f n U ∆φ and σ n P f n P ∆φ, respectively. Following remark 1, a unique parametric increment ∆φ is defined for both the velocity and the pressure mode. To compute ∆φ, an algebraic problem is obtained via the restriction of the high-dimensional equations (12) to the parametric direction I multiplying the momentum and continuity equations by the spatial functions σ n U f n U and σ n P f n P , respectively.
The resulting parametric problem, solved by means of a collocation approach, is where the residuals r n U and r n P of the momentum and continuity equations in the parametric space are given by being a k , k = 0, . . . , 9 a set of coefficients which depend on user-defined data and on previously computed spatial modes as reported in B.

Remark 6.
Contrary to the parametric problem in [45], the second-order term (∆φ) 2 is maintained in equation (17a). Although this term was negligible in laminar simulations, it has been verified numerically that its presence improves the stability of the solution of the parametric step of the Navier-Stokes equations in the turbulent regime.
The alternating direction algorithm in the PGD-NS routine stops when the relevance of the computed increment is negligible with respect to the amplitude of the corresponding mode. Similarly, the global enrichment procedure stops when the contribution of the current mode in the PGD expansion, measured by means of its relative amplitude, is negligible.
Remark 7. As shown in [45], the a priori PGD algorithm devised starting from the separated formulation of problem (12) provides a stable approximation for both velocity and pressure, without the need for tailored pressure corrections required by a posteriori ROMs for incompressible flows [61,62].

Proper generalised decomposition of the turbulence model
The construction of a separated expression ν q t,PGD of the turbulent viscosity using the SA model first requires determining the approximation ν m PGD of the eddy viscosity via a dedicated PGD solver. For this purpose, the integral form of the steady Spalart-Allmaras turbulence model in the high-dimensional space Ω × I is considered, namely for each cell V i , i=1, . . . , N in the computational domain.
The PGD solver for the SA equation is obtained by inserting the separated approximation (7b) of the eddy viscosity into equation (18). Following the rationale utilised for the flow equations, the first term ν 0 PGD in the PGD expansion is selected to enforce the inhomogeneous Dirichlet boundary conditions on Γ in ∪ Γ w . More precisely, Dirichlet data for ν PGD are selected as full order solutions of the SA equation computed using the boundary condition modes of the velocity field. The following modes are determined using a greedy algorithm named PGD-SA: for each new mode σ m ν f m ν ψ m + δ ν m PGD , the prediction σ m ν f m ν ψ m is computed by performing a first iteration of the alternating direction algorithm, whereas the following iterations allow to devise the corresponding correction term σ m ν δ ν m PGD . It is worth recalling that OpenFOAM employs a segragated approach for the simulation of turbulent flows, by solving independently the flow problem and the SA equation. More precisely, the computation of the eddy viscosity is performed using the velocity field U obtained from the RANS equations as input for the solver. The proposed PGD strategy mimics this approach and replaces the velocity field U by its rank-n PGD approximation U n PGD constructed using the procedure described in section 3.  (18), it follows that the unknown increment σ m ν δ ν m PGD can be obtained by solving where the residual R ν is given by Problem (19) can thus be efficiently solved by alternatively restricting it to the tangent manifold associated with the spatial and the parametric coordinates as described in the following sections.

Remark 8.
Similarly to the observation in remark 5, the first, second, fourth, sixth and seventh term in equation (19) mimic the original SA turbulence model (18). Since the remaining integrals are not accounted for in the full order model, and in order to minimise intrusiveness with respect to OpenFOAM core routines, these terms are relaxed during the spatial iteration of the PGD algorithm, by evaluating them using the last computed iteration. This approach is detailed in section 3.4.1.

PGD-SA: the spatial iteration
First, the convective velocity field is replaced by its separated approximation U n PGD (x, µ)= n j=1 σ j U f j U (x)φ j (µ). The SA equation (19) in the high-dimensional space Ω×I is then restricted to the spatial direction multiplying it by ψ m . Moreover, the value of the parametric function is set to ψ m and, from equation (9), the increment σ m ν δ ν m PGD simplifies to ψ m σ m ν ∆f ν . The increment σ m ν ∆f ν thus acts as unknown of the spatial iteration of the PGD procedure for the parametric SA equation. More precisely, a cell-by-cell constant approximation where R m ν denotes the spatial residual of the SA equation and the coefficients β k , k=1, . . . , 7 depend on user-defined data and parametric functions as described in B. Moreover, the production, S m , and destruction, f m w , coefficients are evaluated using the m previously computed modes of ν m PGD and S x PGD ( S m ) and S x PGD (f m w ) denote the corresponding spatial modes of these functions.
As commented in remark 8, three terms appearing in the high-dimensional problem (19) do not have a counterpart in the original SA equation available in OpenFOAM. In order for the discussed PGD-ROM implementation to be non-intrusive with respect to the SA solver natively implemented in OpenFOAM, these terms are treated explicitly via a relaxation approach. Hence, the last three terms on the right-hand side of equation (21) where the parametric function ψ m in the integrals above stems from the restriction of the high-dimensional residuals R ν to the tangent manifold in the spatial direction. By eploiting the separated expression of ν m−1 PGD in terms of its spatial and parametric modes, the residual can be rewritten as where the coefficients β k , k=8, . . . , 12, reported in B, can be precomputed for an efficient implementation of the PGD spatial iteration since they depend soley on user-defined data and parametric functions.

Remark 9.
It is worth emphasising that neither S(x, µ) nor f w (x, µ) are separable exactly via an analytical procedure. It follows that equation (21) cannot be solved in a complete non-intrusive way using OpenFOAM. The resulting implementation of the PGD-SA algorithm in OpenFOAM is thus minimally intrusive as the structure of equation (21) is the same as the original SA equation (6) but a tailored numerical strategy is required to efficiently compute the coefficients and the integral terms involving S and f w . Suitable procedures include a posteriori PGD separation [55,56] and high-dimensional reconstruction of the functions in the space Ω × I followed by an interpolation step in Ω. The latter strategy is employed for the simulations in section 4.

PGD-SA: the parametric iteration
The parametric iteration of the turbulence model is devised by fixing the newly computed spatial mode σ m ν f m ν ←σ m ν [f m ν +∆f ν ] of the eddy viscosity. The corresponding parametric increment ∆ψ is obtained by restricting the high-dimensional equation (19) to the parametric direction I multiplying it by the spatial function σ m ν f m ν . It follows the algebraic equation in which b k , k = 1, . . . , 7 represent the precomputed coefficients reported in B, r m ν is the residual of the SA equation in the parametric space and S µ PGD ( S m ) and S µ PGD (f m w ) denote the parametric modes of the non-separable functions S m and f m w , respectively. It is worth recalling that S m and f m w are evaluated using the m previously computed modes of ν m PGD and require appropriate numerical separation strategies to be computed, see remark 9. Finally, the right-hand side r m ν of equation (24) has the form with the coefficients b k , k = 8, . . . , 12 depending upon data provided by the user and spatial functions, see B.
As for the PGD-NS routine, the alternating direction scheme in the PGD-SA algorithm stops when the computed increment is negligible with respect to the amplitude of the mode. Moreover, the PGD enrichment procedure finishes when the relative amplitude of the current mode is negligible with respect to the previously computed ones.

PGD-ν t : devising a separated turbulent viscosity
The approximation ν m PGD of the eddy viscosity provided by the PGD-SA routine is employed to construct a separated representation ν q t,PGD of the turbulent viscosity, according to the relation (2). This procedure is performed by the PGD-ν t routine, before going back to the PGD-NS algorithm for the RANS equations.
It is worth noticing that the function f v1 in equation (2) is not separable analytically. Hence, as detailed in remark 9, either a numerical PGD separation [55,56] or a highdimensional reconstruction of the function in the space Ω × I to perform interpolation in Ω and I separately, is required. The latter strategy is employed for the simulations in section 4.
For the sake of readability, consider f m v1 obtained using the m computed modes of ν m PGD . The PGD separated expression of this function is given by where S x PGD (f m v1 ) and S µ PGD (f m v1 ) denote the spatial and parametric modes of the function f m v1 , respectively, obtained by means of either of the numerical strategies proposed above.
The separated expression ν q t,PGD reported in equation (7c) can thus be devised in terms of elementary arithmetic operations of separated functions [56], exploiting the separated nature of all the involved quantities. More precisely, it holds where the spatial modes f j t can be computed as the product of the spatial functions S x PGD (f m v1 ) and f m ν , whereas the parametric terms ξ j can be obtained from the product of the parametric modes S µ PGD (f m v1 ) and ψ m . It is worth noting that the product of separated functions leads to a number of terms in the PGD expansion larger than the number of modes of its factors. Nonetheless, such operation can be efficiently performed via a separated multiplication [56] and the result can be compressed to eliminate redundant information [55,56].

Numerical experiments
In this section, numerical simulations of the NASA wall-mounted hump are presented to demonstrate the potential of the proposed methodology. This problem is a quasi-twodimensional NASA benchmark devised to validate turbulence modelling, starting from an experimental set-up. The domain consists of a Glauert-Goldschmied type body mounted on a splitter plate between two endplates, see figure 1(a). Following [63,64], the characteristic length of the problem is set equal to the chord length of the hump c=0.42 m, whereas its maximum height is 0.0537 m and its span is 0.5842 m. Flow separation along the hump is controlled via a suction jet acting through an opening located at 65% of the chord c, as detailed in figure 1(b). In the experimental set-up the opening is connected to a plenum, on the bottom of which suction pumps are installed; for the numerical simulations in the following sections, the plenum is removed and the suction effect is imposed as a boundary condition on the opening via a mass flow rate of 1.518 × 10 −2 kg/s for the jet.
In the analysis of this problem, the quantity of interest is represented by the effect of the suction jet on the flow separation and on the position of the reattachment point. Experimental and numerical studies [65, 66] verified the quasi-two-dimensional nature of the phenomena identifying minor three-dimensional effects located near the endplates. Henceforth, the PGD results will be compared to the full order OpenFOAM approximation, considered as reference solution.  The NASA wall-mounted hump problem being quasi-two-dimensional, in the following sections both the 2D and the 3D cases are studied. A parametric flow control problem is devised by varying the maximum amplitude of the suction jet between 10% and 100% of the module of a peak velocityÛ . In two dimensions, a sinusoidal velocity profile is defined as wherex is the normalised curvilinear abscissa of the jet patch, that isx ∈ [0, 1], and the resulting jet is pointing in the directionŷ orthogonal to the boundary. Similarly, in the 3D case the jet defined on the plane (x,ẑ) and pointing in the orthogonal directionŷ is where the normalised coordinateẑ iŝ It is worth noting that the module of the peak velocityÛ is selected such that the ratio between the mass flow rate of the jet and of the inlet is 1.5 × 10 −3 , reproducing the value in the experimental set-up of the NASA wall-mounted hump with a plenum [63]. In addition, both in (28) and (29), the interval of the parametric variable is defined as I=[0.1, 1].

Two-dimensional NASA wall-mounted hump with parametrised jet
The computational domain for the two-dimensional NASA wall-mounted hump is a channel of height c, extending 6.39c upstream and 5c downstream as displayed in figure 2. The resulting mesh consists of 114, 000 cells. Homogeneous Dirichlet boundary conditions are imposed on both the velocity and the eddy viscosity on the bottom wall and on the hump. A symmetry condition is imposed on the top wall, whereas on the outlet free traction is enforced. At the inlet, a parabolic profile is imposed for both velocity and eddy viscosity. More precisely, the variables range between a null value on the bottom wall and a maximum value at y=0.02 m. For the velocity, the peak value is 34.6 m/s, whereas the free-stream eddy viscosity ν=3ν is selected. The kinematic viscosity is ν=1.55274 × 10 −5 m 2 /s, thus the resulting Reynolds number is approximately Re=936, 000. On the jet patch, the velocity profile (28) withÛ =23.4 m/s is imposed and a homogeneous Neumann condition is considered for the eddy viscosity. It is worth noting that on the hump the mesh is such that y + < 1, whence no wall treatment is required for the turbulent viscosity.
In order to impose the inhomogeneous Dirichlet boundary data, two modes are introduced. In particular, recall that the boundary condition on the jet patch depends upon the parameter µ. For this purpose, two spatial modes, associated with the extreme values of the parametric interval I=[0.1, 1], are computed using simpleFoam with SA turbulence model. The corresponding parametric functions are selected such that to define a linear variation of the Dirichlet data between the extreme values associated with the computed spatial modes. More precisely, the first spatial mode is a full order solution with the jet acting at 100% of the mass flow rate (µ=1) and the second one is associated with the jet acting at 10%, that is µ=0.1. The corresponding parametric modes are given by The tolerance for the enrichment of the flow variables is set to η u =η p =10 −4 , whereas the tolerance for the turbulence model is selected as η ν =10 −2 . The criterion to update the turbulent viscosity is detailed in remark 3, with γ=1. The turbulent pgdFoam algorithm achieves convergence with eight velocity-pressure modes computed using PGS-NS and three corrections by means of PGD-SA and PGD-ν t . Each PGD-SA loop reached the prescribed tolerance within three computed modes. The relative amplitude of the computed modes, as the turbulent viscosity is updated, is reported in figure 3. Following [45], a measure of the relative amplitude accounting for both velocity and pressure modes is employed via the definition It is worth recalling that each time the relative amplitude of the modes drops by one order of magnitude, the separated representation of the turbulent viscosity is updated via the PGD-SA and PGD-ν t routines (see remark 3) and the PGD approximation for velocity and pressure is recomputed using the updated turbulent viscosity.
The importance of updating the PGD approximation of the turbulent viscosity to correctly compute the turbulent velocity and pressure fields is displayed in figure 4(a). Considering the result of the full order simpleFoam with SA turbulence model for µ=0.5 as a reference solution, figure 4(a) compares the relative L 2 (Ω) error of the PGD approximation computed via the PGD-NS, PGD-SA and PGD-ν t strategy described in algorithm 1 with the one obtained by omitting the turbulent viscosity update. Without the PGD-SA and PGD-ν t routines, the error of both velocity and pressure approximations stagnates from the first computed mode and the overall value is one order of magnitude larger than the one achieved by the methodology in algorithm 1. Figure 4(b) reports the relative L 2 (Ω) error of the PGD approximation with turbulent viscosity update for three configurations, that is µ=0.25, µ=0.5 and µ=0.75. The results clearly display that the PGD approximation achieves comparable accuracy throughout the parametric interval I using two boundary condition modes and three computed modes. The following modes only introduce minor corrections to the solution as identified by their corresponding amplitudes, see figure 3.
As mentioned in the problem statement, the quantities of interest in this study are the position of the reattachment point and the effect of the suction jet on the recirculation Qualitative comparisons of the pressure field and the turbulent viscosity for different values of the parameter µ are presented in figure 6 and 7, respectively. Using eight computed modes, the PGD approximation is able to accurately approximate localised variations in the flow pattern, throughout the interval I.
In addition, the accuracy of the PGD-ROM is evaluated by quantitatively comparing µ=0.5 µ=0.75    Similarly to the two-dimensional case, the boundary conditions are enforced using the two parametric modes in (32) and two spatial modes corresponding to the simpleFoam solutions with SA turbulence model for µ=0.1 and µ=1.
The values η u =η p =0.5 × 10 −3 and η ν =10 −2 are considered for the tolerance of the enrichment loops of the flow variables and the turbulent viscosity, respectively. To reduce the overall cost of the PGD-NS, PGD-SA and PGD-ν t procedure, the number of turbulent viscosity updates is reduced by considering a lower initial tolerance in criterion (10), namely γ=2.
Algorithm 1 achieves convergence with four modes computed by the PGS-NS routine and two PGD-SA and PGD-ν t corrections. Each PGD-SA loop reached the prescribed tolerance within two computed modes. The PGD approximation is then compared with the corresponding full order solution provided by simpleFoam with the SA turbulence model: the relative L 2 (Ω) error for µ=0.25, µ=0.5 and µ=0.75 is displayed in figure 10, reporting that the reduced order model is able to provide errors in velocity and pressure below 0.1% and 0.5%, respectively. The resulting PGD-ROM is thus employed to analyse the physical phenomena involved in the turbulent flow over the hump. Figure 11 displays the velocity profile on the hump, computed using the PGD, for different values of the parameter µ. In addition, a qualitative comparison of the streamlines and the recirculation effects computed using the reduced model and the full order OpenFOAM solution are reported in figure 12. The results display that the recirculation effects are reduced when increasing the suction jet and the PGD is able to capture the vortex structure with comparable accuracy with respect to the full order solution.
The capability of the proposed PGD-ROM strategy to treat complex problems of engineering interest is thus confirmed by the following analysis focusing on relevant physical quantities. In figure 13, the top view of the wall shear stress on the bottom wall is reported in the region from the jet patch up to 1.6c downstream, highlighting the effect of the suction jet on flow recirculation. A qualitative comparison between the reduced order and the full order solution confirms the ability of the PGD to accurately reproduce the turbulent flow in the entire range of values I of the parameter. In particular, these conclusions hold true both for the flow variables and for given physical quantities computed starting from them.
Finally, the position of the reattachment point in correspondance of the location of the peak of the jet profile is reported in table 2 for different values of the parameter µ. For each tested configuration, the PGD approximation with four computed modes shows excellent agreement with the full order solver, with relative errors below 0.5%. The results thus display the capability of the PGD-ROM strategy to devise a surrogate model for a quantity of physical interest, robust throughout the parametric domain I, and with an accuracy acceptable for industrial applications.

Conclusion
A PGD strategy to compute parametric solutions of turbulent incompressible flow problems in OpenFOAM has been proposed. The methodology is based on the incompressible Reynolds-averaged Navier-Stokes equations with Spalart-Allmaras turbulence model and mimics the segragated approach implemented in the industrially-validated solver Open-FOAM to devise a minimally intrusive PGD-ROM for convection-dominated flow problems of industrial interest. First, the velocity and pressure modes are computed using the non-intrusive PGD strategy PGD-NS developed in [45] using a seed value for the turbulent viscosity. The PGD approximation of the velocity is then used to improve the turbulent viscosity representation via the minimally intrusive PGD-SA and PGD-ν t routines. Finally, the resulting separated turbulent viscosity is utilised to recompute the PGD expansions of velocity and pressure. The importance of an accurate approximation of the turbulent viscosity has been verified by comparing the solution of the above algorithm with the one computed without solving the SA equation: the latter solution quickly stagnates providing errors of one order of magnitude larger than the proposed methodology.
The developed strategy has been validated in two and three spatial dimensions using a benchmark problem of turbulent external flow, the NASA wall-mounted hump, with Re=93, 600 and Re=936, 000. A flow control problem of industrial interest has been devised by introducing a suction jet on the hump to reduce the recirculation effects. The proposed PGD-based reduced order model has proved to be able to compute a reduced basis with no a priori knowledge of the solution, for convection-dominated viscous incompressible flows achieving both qualitative and quantitative agreement with the full order solution computed via simpleFoam with SA turbulence model, throughout the interval of the parametric variable. More precisely, the reduced model provided accurate approximations of the velocity and pressure fields, with relative L 2 errors below 0.1% and 1%, respectively. In addition, it proved to be able to capture localised flow features and estimate quantities of engineering interest with errors below 0.5%. The reported results thus highlight the robustness of the proposed PGD methodology in presence of turbulent phenomena and its capability to devise accurate approximations of the physical variables involved in the parametric problem (i.e. velocity and pressure), the variables modelling turbulent effects (i.e. eddy and turbulent viscosity), as well as relevant physical quantities, including the position of the reattachment point, the skin friction coefficient and the pressure coefficient.
Finally, it is worth noticing that the minimally intrusive nature of the proposed method represents a promising starting point for the construction of PGD-ROM strategies based on CFD software validated by the industry, beyond incompressible flows. This may be of interest e.g. for parametric compressible flow problems with applications to the aerospace industry. In this context, future investigations will have to deal with the additional difficulty of the treatment of shock waves whose position and intensity may depend upon the value of the parameters of the problem under analysis. A Classical and predictor-corrector PGD algorithms In this appendix, a comparison of two PGD strategies, the classical one [34] and the one based on a predictor-corrector approach [45], is presented using an abstract variational framework with the goal of highlighting main differences, advantages and disadvantages of each solution.
Consider a, possibly nonlinear, partial differential equation (PDE) whose variational form is: seek v(x) ∈ V such that A(w, v(x)) = (w), ∀w ∈ V, where v(x) is the unknwon solution belonging to an appropriately defined functional space V, w is a test function and A and account for the differential operator and the independent term of the problem, respectively.
Given the set of parameters µ ∈ I ⊂ R M , the solution of the corresponding parametric PDE is given by the function v(x, µ) ∈ V µ := V ⊗ L 2 (I 1 ) ⊗ · · · ⊗ L 2 (I M ) satifying A µ (w, v(x, µ)) = L µ (w), ∀w ∈ V µ , First, the classical PGD algorithm for the computation of the separated solution of equation (35) is recalled. The PGD approximation of the high-dimensional unknown function v(x, µ) is thus defined as where f n v (x) and ϕ n (µ) represent the n-th spatial and parametric modes, respectively and σ n v is its corresponding amplitude. It is worth noticing that the previously introduced modes are such that f n v = ϕ n =1, whereas the computed modes before the normalisation procedure are denoted byf v n (x) andφ n (µ).
Algorithm 2 reports the flowchart of the classical PGD algorithm described in [34]. A greedy strategy is employed to compute the terms in the rank-n PGD approximation. More precisely, the algorithm assumes that the approximation v n−1 In addition, the computed predictions before the normalisation procedure are denoted bỹ f v n (x) andφ n (µ) and the corresponding corrections are represented by ∆f v (x) and ∆φ(µ).
The predictor-corrector PGD strategy (see algorithm 3) thus splits the computation of each new mode in two stages: first, the prediction step solves one parametric problem (Algorithm 3 -Line 5) and one spatial problem (Algorithm 3 -Line 8) to computeφ n (µ) andf v n (x), respectively; second, the correction stage performs multiple iterations of the alternating direction method to determine the best corrections ∆φ(µ) (Algorithm 3 -Line 11) and ∆f v (x) (Algorithm 3 -Line 14), given the base solution represented by the predicted modesφ n (µ) andf v n (x).

Algorithm 3 Predictor-corrector PGD algorithm
Require: Stopping criterion η v for the PGD enrichment. Stopping criterion AD stopCrit for the alternating direction scheme. 1: Compute boundary condition modes. Set the value of the spatial modef v n .

7:
Set the value of the normalised parametric prediction ϕ n .

8:
Compute the spatial predictionf v n by solving A µ (w,f v n ϕ n ) = L µ (w) − A µ (w, v n−1 PGD ), projected along the spatial direction.

9:
while AD stopCrit not fulfilled do 10: Set the value of the spatial modef v n .

13:
Set the value of the normalised parametric mode ϕ n . 14: Compute the spatial correction ∆f v by solving A µ (w, ∆f v ϕ n ) = L µ (w) − A µ (w, v n−1 PGD ) − A µ (w,f v n ϕ n ), projected along the spatial direction.

15:
Compute the amplitude σ n v = f v n + ∆f v .

16:
Normalise the spatial mode f n v =f Update the mode counter n ← n + 1.

19: end while
It is worth noticing that the predictor-corrector strategy in algorithm 3 requires the solution of one extra set of parametric and spatial problems per mode with respect to the classical PGD procedure in algorithm 2. Nonetheless, in presence of nonlinear problems like the RANS and the SA equations studied in this work, the predictor-corrector strategy provides increased accuracy in the approximation of the high-dimensional nonlinear differential operator A µ , by successively correcting the base solution identified by the predictions. This allows to reduce the number of iterations required by the alternating direction method to converge and the overall cost of the PGD algorithm. On the contrary, numerical experiments have shown that the computational gain of the predictor-corrector PGD is negligible when linear problems are treated. Finally, the predictor-corrector PGD algorithm also provides a natural stopping criterion for the alternating direction scheme, when the norm of the computed corrections is small with respect to the amplitude of the mode under analysis.