Hybridisable discontinuous Galerkin solution of geometrically parametrised Stokes flows

This paper proposes a novel computational framework for the solution of geometrically parametrised flow problems governed by the Stokes equation. The proposed method uses a high-order hybridisable discontinuous Galerkin formulation and the proper generalised decomposition rationale to construct an off-line solution for a given set of geometric parameters. The generalised solution contains the information for all the geometric parameters in a user-defined range and it can be used to compute sensitivities. The proposed approach circumvents many of the weaknesses of other approaches based on the proper generalised decomposition for computing generalised solutions of geometrically parametrised problems. Four numerical examples show the optimal approximation properties of the proposed method and demonstrate its applicability in two and three dimensions.


Introduction
Reduced order models (ROMs) have become commonplace in many areas of computational sciences and engineering. 40 Some popular ROMs used to reduce the complexity of high dimensional problems include the reduced basis method, 42 the proper orthogonal decomposition (POD) 6,7,29 and the proper generalised decomposition (PGD). [10][11][12] One of the main attractive properties of the PGD is its ability to build reduced basis without prior knowledge of the solution. [10][11][12] However, the intrusive implementation and the difficulty in handling geometrically parametrised problems has often been considered a difficulty when considering its application to complex problems. In recent years, there have been an increase in non-intrusive implementations of the PGD. 25,49,53 In terms of geometrically parametrised problems, early work focused on solutions tailored to specific problems 8,13,26,28 or strategies only applicable in a context of low order approximations. 3,52 More recently, a general approach to deal with geometrically parametrised problems in a CAD environment was proposed. 47 The PGD strategy presented in 47 used a classical finite element (FE) discretisation of Stokes flow problems, leading to the need to use the so-called high-order PGD projection 30 to separate some terms of the weak formulation.
In this work a PGD strategy is proposed in the framework of the hybridisable discontinuous Galerkin (HDG) method. 14,16,22,23,45 The use of a mixed formulation is shown to be beneficial as all the terms of the weak formulation can be written in a separated form, as required by the PGD, without invoking to the memory intensive high-order PGD projection. The use of the HDG method for the spatial discretisation also guarantees that equal order of approximation can be used for all the variables circumventing the so-called Ladyzhenskaya-Babuška-Brezzi (LBB) condition. This is of special importance in this work, where geometrically parametrised domains are considered with curved boundaries. The use of the same degree of approximation for all the variables means that standard isoparametric elements can be used. In contrast, the work in, 47 employing standard FEs, required the use of sub-parametric or super-parametric formulations in the presence of curved boundaries due to the different degree of approximation used for the velocity and pressure, as required to satisfy the LBB condition. Furthermore, the proposed HDG-PGD approach facilitates the imposition of the Dirichlet boundary conditions as in the HDG context all boundary conditions are weakly imposed.
The formulation is presented using Stokes flows as the model problem. However, it is worth mentioning that there has been a substantial effort in developing HDG methods for a variety of problems in different areas of science and engineering 24, 27, 35-37, 39, 44, 48 and therefore, the proposed approach can be easily extended to a wide range of problems. It is also worth noting that the integration within a CAD environment proposed in 47 is also feasible given the recent development of a coupled HDG-NEFEM formulation for fluid 46 and solid mechanics. 43 The structure of the remainder of the paper is as follows. Section 2 presents the Stokes flow problem on a geometrically parametrised domain and the corresponding multidimensional parametric problem. The HDG formulation for the multi-dimensional parametric Stokes problem is described in section 3. The proposed PGD rationale is described in detail in section 4. Section 5 presents a series of numerical examples involving Stokes flow problems in two and three dimensions. Finally, section 6 presents the conclusions of the work that has been presented.
2 Problem statement 2.1 The Stokes problem on a parametrised domain Let us consider a parametrised domain Ω µ ⊂ R nsd , where n sd is the number of spatial dimensions and µ ∈ I ⊂ R npa is a set of geometric parameters that controls the boundary representation of the domain, with n pa being the number of geometric parameters. It is worth noting that the set of geometric parameters can be written as I := I 1 ×I 2 ×· · ·×I npa with µ j ∈ I j for j = 1, . . . , n pa .
For any set of parameters µ, the goal is to find the parametric velocity, u(x µ ), and pressure, p(x µ ), fields that satisfy the Stokes problem given by where ν > 0 is the kinematic viscosity, s is the volumetric source and n µ is the outward unit normal vector to ∂Ω µ . The boundary of the domain, ∂Ω µ , is partitioned into the non-overlapping Dirichlet, Γ µ D , Neumann, Γ µ N , and slip, Γ µ S , boundaries such that On the Dirichlet boundary the velocity is given by u D . On the Neumann boundary the pseudo-traction is given by g N . Finally, on the slip boundary, the matrices D µ and E µ are given by D µ = [n µ , 0 nsd×(nsd− 1) ] and E µ = [0, t µ 1 , ..., t µ nsd−1 ], as detailed in. 23 The tangential vectors t µ k , for k = 1, . . . n sd − 1 are such that {n µ , t µ 1 , ..., t µ nsd−1 } form an orthonormal system of vectors.
The free divergence condition in equation (1) induces the compatibility condition where ·, · S denotes the standard L 2 scalar product in any domain S ⊂ ∂Ω µ .
In addition, it is worth noting that, if Γ µ N = ∅, an additional constraint to avoid the indeterminacy of the pressure is required. One common option 15,18,21,32 that is considered here, consists of imposing the mean pressure on the boundary of the domain, namely

The multi-dimensional parametric Stokes problem
The classical strategy to solve the parametric Stokes problem is to solve equation (1) for every set of parameters µ ∈ I. However, this strategy is not well suited when fast queries are required.
Reduced order models have demonstrated to be a viable alternative to compute multidimensional parametric solutions in an offline phase. Once the offline solution is available, the computation of the solution for a given set of parameters has a very small computational cost, being very well suited for applications where fast queries are required.
The multi-dimensional parametric problem arises from interpreting µ as additional parametric coordinates, rather than parameters of the problem. In the context of the Stokes problem considered here, the strategy is to consider the velocity and pressure fields as functions in a multidimensional space, namely u(x µ , µ) and p(x µ , µ). The multidimensional parametric Stokes problem can be written as For the multi-dimensional problem, the compatibility condition induced by the free divergence condition can be written as and the additional constraint to avoid the indeterminacy of the pressure, required when

Hybridisable discontinuous Galerkin formulation
Let us consider a subdivision of the domain Ω µ in n el disjoint subdomains Ω µ e such that The interior boundaries of the subdomains define the so-called mesh skeleton or internal interface Γ µ as A partition of the parametric domains I j , for j = 1, . . . , n pa , in n j el disjoint subdomains I j e such that is also considered.
This section briefly presents the HDG formulation for the multi-dimensional parametric Stokes problem. The presentation is based on previous work on HDG methods found in. 15,17,18,21

Mixed formulation
Introducing the so-called mixed variable L = −ν∇ µ u, the Stokes problem can be written as a first-order system of equations in the broken computational domain, namely in Ω µ e × I, and for e = 1, . . . , n el , ∇ µ · L e + p e I nsd = s in Ω µ e × I, and for e = 1, . . . , n el , ∇ µ · u e = 0 in Ω µ e × I, and for e = 1, . . . , n el , where the last two equations, known as transmission conditions, impose the continuity of the velocity and the normal flux on the mesh skeleton. Following, 31 the jump operator · is defined as the sum from the left, Ω l , and right, Ω r , elements of a given portion of the interface Γ µ × I, that is

Strong form of the local and global problems
The HDG method solves the mixed problem of equation (10) in two steps. First, the so-called local problems are considered in Ω µ e × I, and for e = 1, . . . , n el , ∇ µ · L e + p e I nsd = s in Ω µ e × I, and for e = 1, . . . , n el , ∇ µ · u e = 0 in Ω µ e × I, and for e = 1, . . . , n el , = ρ e , for e = 1, . . . , n el , whereû is the so-called hybrid variable, which is an independent variable representing the trace of the solution on the element faces, and ρ e is the mean value of the pressure on the boundary ∂Ω e . It is worth noting that the local problem is a pure Dirichlet problem and therefore, the last condition in equation (12) is introduced to ensure the uniqueness of the pressure. The local problems can be solved independently, element by element, to write L e , u e and p e in terms ofû and ρ e along the interface Γ µ ∪ Γ µ N ∪ Γ µ S . Second, the so-called global problem is defined to impose the continuity of the normal flux on the inter-element faces and the Neumann and slip boundary conditions, namely It is worth noting that, due to the unique definition of the hybrid variable on each face and the Dirichlet boundary condition in the local problems, there is no need to enforce the continuity of the solution in the global problem.
The constraint of equation (5), induced by the incompressibility condition, is also considered in the global problem and written in terms of the hybrid variable as

Weak form of the local and global problems
The following discrete functional spaces are introduced: where P k (Ω µ e ), P k (Γ µ i ) and P k (I j e ) stand for the spaces of polynomial functions of complete degree at most k in Ω µ e , on Γ µ i and in I j e respectively. The weak form of the local problems, for e = 1, . . . , n el , reads: given where the multi-dimensional bilinear and linear forms of the local problem are given by and respectively, where (·, ·) D denotes the standard L 2 scalar product in a generic subdomain D and τ µ is the stabilisation tensor, whose selection has an important influence on the accuracy, stability and convergence properties of the resulting HDG method. 16,[32][33][34] The choice of the stabilisation tensor for geometrically parametrised problems will be discussed in the next section.
Similarly, the weak form of the global problem is: for allv ∈ V h µ , where the multi-dimensional bilinear and linear forms of the global problem are given by and respectively.

The proper generalised decomposition strategy
The solution of the parametric problem of dimension n sd + n pa , presented in the previous section, with the standard HDG approach is usually not affordable, even for a relatively small number of parameters. To circumvent the curse of dimensionality, this section proposes the use of the PGD framework. As it will be shown in this section, the use of an HDG formulation has important advantages compared to other formulations such as standard finite elements. 47 To simplify the presentation, the subindex e and the superindex h used in the previous section to specify the element and the discrete approximations will be omitted here, unless they are needed to follow the development.

Separated spatial mapping to obtain generalised solutions
As discussed in detail in, 38,41,47 the solution of the parametric problem described in section 3 requires that the bilinear and linear forms in the weak form can be expressed, or well approximated, by a sum of products of parametric functions and operators that are parameter-independent. To enforce the affine parameter dependence, the integrals appearing in the weak form must involve domains that are not dependent upon the parameters. Following the work of, 3, 47, 52 a mapping between a parameter-independent reference domain, Ω, and the geometrically parametrised domain is considered, namely The coordinates of the reference, or undeformed, domain are denoted by x whereas the coordinates of the parametric, or deformed, domain are denoted by x µ . To ensure the affine parameter dependence, the mapping is assumed to be given in separated form as Remark 1. To simplify the presentation here, it is assumed that the separated representation of the mapping is given analytically. As mentioned earlier, a general strategy to construct a separable mapping was described in 47 using an exact boundary description of the computational domain by means of NURBS.
The separated representation of the mapping leads to the following separated representation of its Jacobian In addition, the separated description of the mapping and its Jacobian can be used to obtain a separated expression of the determinant and the adjoint of the Jacobian using the Leibniz formula and the Leverrier's algorithm as explained in detail in. 47 The separated expression of the determinant of the Jacobian and its adjoint are written in compact form as and respectively.

Affine parameter dependence of the HDG bilinear and linear forms
Introducing the mapping M µ of equation (21) into the weak form of the local and global problems, it is possible to write the integrals over the reference domain, Ω, and its boundary, ∂Ω, not dependent on the parameters µ. The bilinear and linear forms for the local problems can be written as and respectively, where the adjoint operator is defined as adj(A) = det(A) A −1 and the stabilisation parameter in the deformed domain is chosen as The scaling factor adj(J µ )n in equation (28) accounts for the increased or decreased area of the deformed face, ∂Ω µ e , with respect to the reference one, ∂Ω e . This definition, inspired by the expression of the penalty coefficient in classical interior penalty DG methods, 4 ensures that the larger the deformation of the face, the smaller the value of τ µ is. This ensures that a weaker continuity is imposed for large deformations and it is justified by the expected loss of accuracy in the hybrid variable when the mapping introduces a large deformation.
Following previous work on HDG methods for Stokes problems, 21 the stabilisation parameter in the reference domain is selected as τ = (τ ν/ )I nsd , where τ is a numerical parameter, selected as τ = 10 in this work, and is a characteristic length of the domain. Remark 2. As mentioned above, it holds that adj(J µ )n = |∂Ω µ e |/|∂Ω e |. Hence, no parametric dependence appears in the arguments of the bilinear form A ρp .
Analogously, the bilinear and linear forms for the global problem can be written as and Lû respectively.

Remark 3.
The derivation of the terms on the slip boundary in (29) follows from the (28). The slip boundary condition is used here to enforce a symmetry condition and therefore, it is assumed that the orientation of the vectors {n µ , t µ 1 , ..., t µ nsd−1 } is preserved by the mapping M µ . It is worth noting that this does not imply that Γ µ S = Γ S as it will be shown with numerical examples.

Remark 4.
As classical in the context of shape optimisation, 1 in (30) it is assumed that Neumann boundaries, where a traction (or pseudo-traction) is imposed, are fixed, that is, Γ µ N = Γ N . On the contrary, deformable Neumann boundaries, also known as free boundaries, are traction-free, whence g N is null.

Separated representation of the data
As usual in a PGD context, the data is assumed to be given in separated form. For the Stokes problem under consideration, this means that the Dirichlet and Neumann data and the source term can be written as Even if the data is not directly given in this form, it is possible to obtain a good approximation in a separated form, see. 11

Separated representation of the primal, mixed and hybrid variables
Following the predictor-corrector PGD rationale, see, 49 each variable of the HDG formulation, presented in section 3, is written as a rank-m separable approximation, that is ρ ∆ρ m PGD are the corresponding correction terms. Introducing the variation ∆, the correctors are defined as where the least term denotes a high-order variation and it is henceforth neglected.
Each term, or mode, of the PGD approximation is the product of a function that depends upon the spatial coordinates and a function that depends upon the parameters. In addition, the parametric functions are assumed to be the product of functions that depend upon a single parameter, namely As usual in a PGD context, the number of terms is a priori unknown and it is determined using a greedy algorithm. Assuming that m − 1 modes are known, the computation of the m-th mode requires the solution of a nonlinear problem as described in the next section.
Remark 5. This work considers the so-called single-parameter approach, where the parametric function of the m-th mode, ψ m , is the same for all the variables. Other approaches, including a different parametric function for each variable or even the use of vector-valued parametric functions in the approximation of vector fields are discussed in. 20 The tangent manifold for L is characterised by choosing W as variations of F L and ψ, that is for δF L ∈ W h := V h (Ω) nsd×nsd and δψ ∈ L h (I). Similarly, the tangent manifolds for u, p, u and ρ are characterised by choosing

Alternating direction scheme
With the separated structure of the PGD approximations, the weighting functions and the bilinear and linear HDG forms described in the previous sections, it is possible to drastically reduce the complexity of the problem by projecting the high-dimensional problem on the tangent manifold and applying an alternating direction strategy. First, in the so-called spatial iteration, the parametric function of the m-th mode is assumed known and the spatial functions are determined. As it will be shown, this step requires to solve a system of equations with a very similar structure to the non-parametric HDG problem. Second, in the so-called parametric iteration, the parametric function is computed using the spatial functions determined in the first step. This process is repeated until convergence is achieved. It is worth noting that the order of the spatial and parametric iterations can be swapped without affecting the alternating direction algorithm.
Let us assume that we have computed the first m − 1 modes and it is of interest to compute the m-th mode. In the next two sections, the alternating direction strategy to compute the spatial and parametric modes is detailed.

The spatial iteration
In the spatial iteration, it is assumed that the parametric function ψ m and the spatial predictions and σ m ρ f m ρ are known and the goal is to compute the corresponding corrections σ m L ∆F L , σ m u ∆f u , σ m p ∆f p , σ m u ∆fû and σ m ρ ∆f ρ . Taking into account that δψ = 0 when ψ m is known and introducing the expression of the PGD approximations and the weighting functions in the weak form of the HDG local problems, the following weak form of the local problem for the spatial iteration is obtained: The bilinear and linear forms of the local problem are detailed in equation (58), in A, and equation (61), in B, respectively. The constants in equation (37) are given by where the bilinear forms involved in the definitions of these constants are introduced in equation (60), in A.
As mentioned earlier, in remark 5, this work considers the same parametric function for all the variables. It is worth noting that this choice reduces the number of different constants in equation (37).
Similarly, the weak form of the global problem is: The bilinear and linear forms of the global problem are detailed in equation (59), in A, and equation (62), in B, respectively.

The parametric iteration
After computing the spatial corrections following the procedure described in the previous section, the spatial modes are updated, namely where the constant σ m on the left hand side denotes the amplitude of the newly computed m-th mode of the function , e.g. σ m p ← σ m p f m p + σ m p ∆f p . In the parametric iteration, the goal is to compute the parametric correction ∆ψ given the prediction ψ m and the known spatial functions in (40). Following the assumption that such functions are known, it holds that δF L = δf u = δf p = δfû = δf ρ = 0. Introducing the expression of the PGD approximations and the weighting functions in the weak form of the HDG local problems, the following weak form of the local problem for the spatial iteration is obtained: find ∆ψ ∈ L h (I) such that for all δψ ∈ L h (I), Similarly, the weak form of the global problem is: find ∆ψ ∈ L h (I) that satisfies for all δψ ∈ L h (I).
The constants in equations (41) and (42) are defined as The choice of a single parameter approximation implies that we can combine equations (41) and (42) to obtain the following parametric problem: find ∆ψ ∈ L h (I) that satisfies for all δψ ∈ L h (I), where Remark 6. Alternative formulations of the parametric problem may be devised, e.g. by considering only equation (41) or (42). In this work, equation (44) has been considered in the parametric iteration in order to account for the information of both the local and the global HDG problems.
As detailed in equation (34), the parametric iteration involves n pa geometric parameters. To reduce the size of the problem of the parametric iteration, n pa one-dimensional problems are solved sequentially, as commonly done in a PGD framework. Set q ← 1 and initialise the parametric predictor ψ m ←1.

13:
Normalise the spatial predictor σ m u ← σ m u f m u + σ m u ∆fû . 14: Update the spatial predictor σ m u f m u ←σ m u f m u + σ m u ∆fû.

15:
Update the stopping criteria for the mode amplitude εû← σ m u ∆fû /σ m u and the residuals ε r • ← r • .

16:
Increase the counter of the alternating direction iterations q ← q + 1.

19: end while
In the greedy enrichment loop, first a predictor of the spatial mode is computed as the solution of the HDG global and local problems using a guess for the parametric mode (Algorithm 1 -Steps 3-5). Then, the alternating direction scheme computes the corrections of the parametric (Algorithm 1 -Steps 8-10) and spatial mode (Algorithm 1 -Steps 11-14) solving a parametric linear system and the HDG global and local problems, respectively. The nonlinear iterations of the alternating direction scheme stop when the amplitude σ m u ∆fû of the correction is negligible with respect to the amplitude σ m u of the current mode and the residuals of the spatial and parametric problems are below a given tolerance (Algorithm 1 -Steps 7 and 15). The stopping criterion for the greedy enrichment algorithm relies on the relative amplitude σ m u of the current mode being negligible with respect to the first mode σ 1 u (Algorithm 1 -Step 2). Alternative stopping criteria based on normalising the amplitude of the current mode with respect to the cumulative amplitudes of the previous modes have also been considered in the literature, see e.g. 49 Note that for the purpose of normalisation (Algorithm 1 -Step 14), an appropriate norm needs to be defined and the L ∞ norm has been utilised for the simulations in section 5.

Discretisation of the spatial and parametric problems
The discretisation of the local problems of the spatial iteration using an isoparametric formulation with equal interpolation for all the variables, 44-46 leads to a system of equations for each element with the following structure: where F L , F u , F p and Fû denote the nodal values of the unknown spatial functions σ m L ∆F L , σ m u ∆f u , σ m p ∆f p and σ m u ∆fû respectively and the constraint on the mean value F ρ of the pressure on the element boundaries is enforced using the Lagrange multiplier F ζ .
The only difference between the local system obtained in the spatial iteration of the proposed HDG-PGD approach and the local system of a standard HDG method 44,46 lies in the construction of the blocks forming the matrices A and vectors f . As an example, let us consider the matrix A LL . In the proposed HDG-PGD framework, this matrix is defined as whereas in a standard HDG approach, the corresponding matrix is defined as In the above expressions {N I } denotes the set of shape functions used to define the spatial approximation of the mixed variable.
Similarly, the discretisation of the global problem of the spatial iteration leads to a system of equations for the trace of the velocity on the element boundaries and the mean value of the pressure in each element, namely As usual in an HDG context, the local problem of equation (46) is used to express the spatial part of the gradient of the velocity, the velocity and the pressure in terms of the spatial part of the trace of the velocity and the mean pressure. Introducing these expressions into the global problem, leads to the global system where the only unknowns are the spatial parts of the trace of the velocity and the mean pressure.
In a similar fashion, the discretisation of the parametric problem (44) using Lagrange shape functions leads to an algebraic system of equations whose unknowns are the nodal values of the parametric modes.

A remark for a computationally efficient implementation
The evaluation of the right hand sides of the PGD spatial and parametric iterations tends to become computationally expensive when approximations with a large number of modes are considered. Indeed, the number of terms involved in such computation experiences a geometric growth rate during the iterations of the greedy algorithm.
In order to ease the computational burden of the overall algorithm, the number of terms in the modal approximations u m PGD , p m PGD , L m PGD ,û m PGD and ρ m PGD is reduced. It is well known that the terms in the PGD reduced basis are not orthogonal to each other and repeated information may appear. Hence, orthogonal separable approximations featuring m < m modes are constructed via the PGD compression, 19, 30 that is, a least-squares higher-order projection minimising the L 2 norm of the difference between target and test functions, From a practical point of view, the PGD compression is applied during the enrichment strategy described in algorithm 1. A trade-off between the cost of performing the greedy iterations with a larger number of modes and the extra cost required by the PGD compression needs to be achieved. For the simulations in section 5, PGD compression is applied every ten new computed modes for the analytical examples and every five for the microfluidics test cases.

Coaxial Couette flow
The first example considers the well known coaxial Couette flow problem, 9 consisting of the flow confined within two infinite coaxial circular cylinders with radius R in and R out respectively, with R in < R out . The boundary conditions introduce the known angular velocities, Ω in and Ω out , at R in and R out , respectively. The problem has analytical solution, given by the azimuthal component of the velocity as where r is the distance to the axis of the cylinders.
To demonstrate the applicability of the proposed ROM the problem is considered in two dimensions, with Ω µ = {x µ ∈ R 2 | µ 1 ≤ r µ ≤ R out }, with R out = 5 and µ 1 ∈ I = I 1 = [1, 3] and where r µ = (x µ 1 ) 2 + (x µ 2 ) 2 . The reference domain is chosen to be Ω = {x ∈ R 2 | 1 ≤ r ≤ R out } and the mapping between the reference and the geometrically parametrised domains is defined by the general separable expression of equation (21) with the mapping of equation (22) given by where r = x 2 1 + x 2 2 . The Jacobian of the mapping is also written in the general separated  form of equation (23), with For the numerical experiments in this section, four triangular meshes of the reference domain are generated, as shown in Figure 1. The meshes have 128, 512, 2,048 and 8,192 elements respectively.
The proposed HDG-PGD framework is used to obtain the generalised solution of the parametric Stokes problem. The first four normalised modes of the magnitude of the velocity field are displayed in figure 2. The computation was performed using the second mesh shown in figure 1 with a degree of approximation k = 4 for all the variables and with a mesh of 1,000 elements in the parametric dimension with also k = 4. As usual in a the context of ROMs, the first modes capture the most relevant and global features of the solution whereas the features captured for the next modes only introduce localised features. Figure 3 shows the first eight normalised parametric modes computed. It can be observed that the first three modes are smooth, whereas the next modes, that have a less relevant contribution to the generalised solution, show a more oscillatory character.  To quantify the importance of the modes on the generalised solution, figure 4 shows the relative amplitudes of the modes with respect to the amplitude of the first mode for all the variables. It can be clearly observed that the fourth mode has an amplitude that is already more than 100 times smaller than the amplitude of the first mode. After computing only nine modes the relative amplitude is already of the order of 10 −6 . It is worth noting that in practice it is not required to add modes with such a lower relative amplitude with respect to the first mode, but in this first example nine modes are computed to show the rapid decrease in their amplitudes.
Once the generalised solution is computed, it is of interest to quantify its accuracy. Figure 5 shows the absolute value of the error of the velocity magnitude using as the number of modes is increased for three relevant configurations corresponding to the parameter µ 1 = 1, µ 1 = 2 and µ 1 = 3. The results show that with only one PGD mode an absolute error below 10 −1 is already obtained for all three configurations, with more accurate results for the case with µ 1 = 2. With two PGD modes the error drops substantially, being less than 7 × 10 −3 in all cases, and with only PGD modes the error is below 2 × 10 −4 for the three configurations considered.
To further illustrate the accuracy of the proposed HDG-PGD approach, the relative error in the L 2 (Ω × I) norm, defined as is studied and compared to the error of the full order HDG approach. Figure 6 shows the evolution of ε PGD , for all the variables, as the number of PGD modes is increased, for different meshes using a quadratic degree of approximation. The discontinuous lines in Figure 6 show the relative error of the full order HDG method, measured in the L 2 (Ω × I) norm. It is worth noting that the computation of the error for the full order approach requires the computation of a large number of solutions. More precisely, the number of HDG solutions required is equal to the number of elements in the parametric space multiplied by number of integrations points in each element.
The results show that the error of the proposed ROM converges monotonically to the error of the full order approach with as the number of modes is increased. In all cases the number of PGD modes required to reach the maximum accuracy on a given mesh is lower than six. Furthermore, the results in figure 6 illustrate the increased level of accuracy obtained as the spatial and parametric discretisations are refined. Analogous results, not reported here for brevity, are obtained for lower and higher orders of approximation.
Next, the optimal approximation properties of the proposed HDG-PDG method are studied by performing a mesh convergent study. Figure 7 shows the evolution of the relative error in the L 2 (Ω × I) norm as a function of the characteristic element size, h, for different orders of approximation and for all the variables of the HDG formulation. The optimal rate of convergence, equal to h k+1 , is approximately observed for all the variables. In each case, the minimum number of PGD modes required to achieve the accuracy of the full order method is selected, as previously discussed when presenting the results of figure 6.
Finally, it is worth mentioning the differences between the proposed HDG-PGD approach presented here and the recently proposed PGD approach for geometrically parametrised domains in 47 using standard finite elements for the spatial discretisation. First, the current approach does not require the higher order PGD projection to separate the inverse of the determinant of the Jacobian, given the first-order character of the problem solved with HDG. Second, the current approach enables the use of the same degree of approximation for velocity and pressure, contrary to the standard FE approach where specific choices are required to satisfy the LBB condition. In the context of geometrically parametrised domains with curved boundaries this implies that the current approach enables the use of isoparametric elements whereas super-parametric or sub-parametric elements are required in the FE context. Third, the weak imposition of the Dirichlet boundary conditions, as usually done in a DG context, facilitates the construction of the generalised solution without the need for specific choices for the modes that satisfy the Dirichlet boundary conditions, as required by approaches. Finally, the results in figure 7 can be compared to the results in. 47

Axisymmetric Stokes flow past a sphere
The second example considers the Stokes flow past a sphere, a typical test case for axisymmetric Stokes flow solvers. The domain of interest is selected as the region confined by two concentric spheres with radius R in and R out respectively, with R in < R out . This problem also has analytical solution, given, in polar coordinates, by the following velocity and pressure fields where v ∞ and p ∞ are the magnitude of the velocity and the pressure of the undisturbed flow, far away from the obstacle. A typical quantity of interest in this problem is the drag force, whose exact value is given by F D = 6πνv ∞ R in Similar to the previous example, the geometric parameter considered here is the radius of the inner sphere. The parametric domain considers the axial symmetry of the problem is defined as Ω µ = {x µ ∈ R 2 | x µ 2 ≥ 0 and µ 1 ≤ r µ ≤ R out }, with R out = 5 and µ 1 ∈ I = I 1 = [1,3]. The reference domain is chosen to be Ω = {x ∈ R 2 | x 2 ≥ 0 and 1 ≤ r ≤ R out }. The mapping between the reference and the geometrically parametrised domains is exactly the same mapping utilised in the previous example, given by the two terms in equation (52).
A no-slip boundary condition is imposed on the inner sphere, a Dirichlet boundary condition corresponding to the exact solution on the outer boundary and axial symmetry is imposed on the rest of the boundary. The axial symmetry is imposed by selecting α = β = 0 in the matrices D µ and E µ in equation (1). As mentioned earlier, in Remark 3,   the portion of the boundary where the axial symmetry is imposed depends on the geometric parameter, but the normal and tangent to the boundary are independent on the geometric changes. Therefore, the matrices D and E do not depend upon the geometric parameters.
The proposed ROM is used to obtain the generalised solution of the parametric axisymmetric Stokes problem. The first four normalised modes of the magnitude of the velocity field and the pressure are shown in figures 8 and 9. The computation was performed using the second mesh with a degree of approximation k = 4 for all the variables and with a mesh of 1,000 elements in the parametric dimension with also k = 4. Figure 10 shows the first eight normalised parametric modes computed. It is worth noting that despite the different Figure 11: Axisymmetric flow past a sphere: Convergence of the mode amplitudes.
nature of the flow and the axisymmetric boundary condition, the parametric modes have a similar behaviour when compared to the modes obtained in the previous example. This is mainly attributed to the geometric parameter describing an analogous variation of the computational domain.
As in the previous example, the evolution of the relative amplitude of the modes is shown in 11. The rapid decrease shows that it is possible to compute a generalised solution to this problem with a very small number of modes. With eight modes the relative amplitude is already below 10 −5 .
Next, the optimal approximation properties of the proposed HDG-PGD method are studied by performing a mesh convergent study. Figure 12 shows the evolution of the relative error in the L 2 (Ω × I) norm as a function of the characteristic element size, h, for different orders of approximation and for all the variables of the HDG formulation. The optimal rate of convergence, equal to h k+1 , is approximately observed for all the variables.
Finally, the accuracy of the HDG-PGD approach on the drag force is studied for three different configurations corresponding to µ 1 = 1, µ 1 = 2 and µ 1 = 3. Figure 13 shows evolution of the error in the drag force as the number of of degrees of freedom is increased for the three different geometric configurations and for different orders of approximation. The number of degrees of freedom refers to the size of the HDG global problem as this is the most time consuming part of the spatial iteration.
The results show the variation of the drag force induced by the variation of the geometric parameter and how the generalised solution produces accurate results for any value of the geometric parameter. In all cases, convergence to exact value is observed, and the superiority of using high-order approximations is clearly appreciated. For the first configuration, the results in figure 13(a) show that with a linear approximation requires the solution of a global problem with 24,832 degrees of freedom to obtain relative error in the drag force of 0.0181. In contrast , using a quartic approximation, the error in the first The results also show that for higher values of the geometric parameter the solution is slightly more difficult to capture and the number of degrees of freedom required is slightly higher. In fact, the advantages of high-order approximations are more noticeable for the case of µ 1 = 3.

Axisymmetric Stokes flow around two micro-swimmers
The next example considers the Stokes flow around the so-called push-me-push-you microswimmer, proposed in. 5 This swimmer consists of two spherical bladders that have the ability to change their mutual distance and individual volume, whilst maintaining the total volume of the two spheres. The swimmer is placed in a cylindrical channel of length L and diameter D.
Two geometric parameters are considered in this example. The first one, µ 1 ∈ I 1 = [−1, 1], controls the radius of the two spheres in such a way that the total volume of the two spheres is maintained. The second parameter, µ 2 ∈ I 2 = [−3, 2], controls the distance between the centre of the two spheres. The value of µ 1 = −1 corresponds to the configuration where the radius of the first sphere is R 1 = 0.3096 and the radius of the second sphere is R 2 = 0.116, whereas the value of µ = 1 corresponds to the opposite situation, with R 1 = 0.116 and R 2 = 0.3096. The value of µ 2 = −3 corresponds to the case where the distance between the spheres is maximum, with the centres of the spheres placed at (−3, 0) and (3, 0) respectively. The value of µ 2 = 2 corresponds to the case where the distance between the spheres is minimum, with the centres of the spheres placed at (−0.5, 0) and (0.5, 0) respectively.
Using the axial symmetry of the problem, the reference domain is chosen as Ω = where L = 6, H = 2, x 0 = (1.5, 0) and R ref = 0.116. Figure 14 shows the triangular mesh of the reference domain used for this numerical example. The mesh has 1,426 elements, leading to a system in the HDG global problem of 22,260 equations for a degree of approximation k = 4.
On the left part of the boundary a Dirichlet boundary condition, corresponding to a horizontal velocity of magnitude one, is imposed. On the right part of the boundary a   homogeneous Neumann boundary condition is imposed. On the surface of the two spheres a no-slip boundary condition is enforced and on the rest of the boundary a slip boundary condition is imposed.
The geometric mapping used in this example is detailed in C.
The first four spatial modes for the velocity and pressure computed with the proposed HDG-PGD are shown in figures 15 and 16. The computation was performed using the mesh of figure 14 with a degree of approximation k = 4 for all the variables and with a mesh of 10,000 elements in each parametric dimension with also k = 4. It is worth noting that the cost of the one-dimensional parametric problems is negligible when compared to the cost of the spatial iteration. Therefore, a large number of elements is used in the parametric dimension to ensure that the variation induced by the geometric parameters are captured with no a priori knowledge of the solution.  figure 17(a) the first, third, fifth and six parametric modes have a normalised value near one for the whole range of values of µ 1 . A similar behaviour is observed for the second parameter µ 2 . In addition, the second parameter, corresponding to the distance between the spheres it can be observed that many of the modes have a much more relevant influence near µ 2 = 2. This is expected as this configuration corresponds to the case where the distance between the spheres is minimum and therefore induces an important variation in the flow field because the first sphere will influence the flow that is reaching the second sphere.
The evolution of the relative amplitude of the modes is displayed in figure 18. The results show that with 24 modes all the relative amplitude of the hybrid variable, used to check convergence, is below 10 −3 . A slower decrease of the relative amplitudes when compared with the previous examples can be observed. This is attributed to two factors. First, this problem considers two geometric parameters and, second, the range of variation of the distance is relatively high when compared to the minimum radius of the spheres.
To illustrate the variation in the geometry induced by the parameters as well as the different flow features that are induced by the geometric changes, figure 19 shows the magnitude of the velocity and the pressure fields in the three dimensional domain for     To analyse the accuracy of the proposed approach, figure 20 compares the drag force on the two spheres as a function of the µ 2 , controlling the distance between the spheres, and for three different configurations of the µ 1 , controlling the radius of both spheres. The results obtained with the HDG-PGD approach are compared to the results of the standard HDG method on a reference mesh. Both solutions show an excellent agreement in all cases.
Finally, to stress the potential of the proposed approach, figure 21 shows the drag force on the two spheres and the total drag as a function of both geometric parameters. This figure shows that generalised solution computed with the HDG-PGD approach can be used to rapidly explore the whole space of parameters and used to find optimal strokes, of interest in many applications. 2

Stokes flow around a sphere in a corrugated channel
The last example, inspired from the studies in, 50, 51 considers the flow past a sphere placed in a corrugated channel. The corrugated channel has a height of 1µm and the undulatory profile is defined by the expression where L = 12.5µm, f ω = 2µm and the value of f n controls the oscillation of the boundary. A sphere of radius R, centred at the origin, is placed inside the corrugated channel. To demonstrate the applicability and potential of the proposed methodology in three dimensions, two geometric parameters are considered. The first parameter µ 1 ∈ [−1, 1] is used to control the radius of the sphere, defined as R(µ 1 ) = (µ 1 + 2)/10. The second parameter µ 2 ∈ [0, 2] controls the amplitude of the corrugated channel, given by f n = 1/2 + µ 2 . The geometry of the reference domain, corresponding to µ 1 = µ 2 = 0, is shown in Figure 22(a). Exploiting the symmetry of the problem, a mesh of a quarter of the domain is considered, with 2,191 tetrahedral elements, as depicted in Figure 22(b).
The geometric mapping used in this example is detailed in D.  example, 12 modes are required to ensure the relative amplitude of the hybrid variable, used to check convergence, is below 10 −3 . Figure 27 shows the magnitude of the velocity and the pressure fields in the channel for three different configurations. The results illustrate the variation in the velocity and pressure fields as the amplitude of the channel and the radius of the sphere is increased.
To assess the accuracy of the computed generalised solution computed with the proposed approach, a reference solution is computed for the three configurations displayed in Figure 27. The reference solutions are computed on a much finer mesh with a standard  HDG solver. As a quantity of interest, the drag on the sphere is measured. Figure 28 shows the evolution of the error of the drag force as the number of PGD modes is increased. To further analyse the accuracy of the computed generalised solution, the error of an HDG solution, computed in each configuration using the same spatial resolution as the one used in the HDG-PGD formulation is considered. The results show that the error of the HDG-PGD approach tends to the error of the HDG solution computed for each configuration, showing the ability of the proposed approach to accurately capture the solution for different geometric configurations.
As mentioned in the previous example, the proposed approach provides a generalised solution that can be used to perform fast queries of different quantities of interest. To illustrate the potential of the developed HDG-PGD approach, Figure 29 shows the drag force on the sphere and the pressure drop, measured as the difference between the pressure at the inlet and outlet, as a function of the geometric parameters µ 1 and µ 2 . The results show that the drag force is not sensitive to the variation of the amplitude of the channel oscillation but very dependent on the radius of the sphere. In contrast, the pressure drop shows a dependency on both geometric parameters.

Concluding remarks
A reduced order model approach based on the PGD and the HDG methods is being presented for the solution of geometrically parametrised Stokes flow problems. The mixed formulation characteristic of HDG methods is beneficial in the PGD context as all the terms of the weak formulation can be written in a separated form, without using to the memory intensive high-order PGD projection. The use of the HDG formulation also enables the use of equal order of approximation for all the variables circumventing the LBB condition. This is advantageous in the context of geometrically parametrised problems in complex domains as it enables the use of standard isoparametric formulations. In addition, the use of a DG formulation implies that no special treatment of the Dirichlet boundary conditions is required.
The optimal approximation properties of the proposed approach have been validated numerically using two and three dimensional test cases. In addition, the ability of the proposed approach to compute generalised solutions involving geometric parameters has been illustrated for problems relevant to the microfluidics community. The examples involve geometric parameters that involve substantial changes of the geometry and induce important changes in the flow features and the relevant quantities of interest.

A Bilinear forms of the HDG-PGD weak formulation
The bilinear forms introduced in the spatial iteration are given by for the HDG local problems and by for the HDG global problems.
In addition, the following bilinear forms are introduced in the parametric iteration

B Linear forms of the HDG-PGD weak formulation
The linear forms introduced in the spatial and parametric iterations are given by for the HDG local problems and by for the HDG global problems.

C Geometric mapping for the channel with two microswimmers
The mapping used in the example involving the flow around two microswimmers is designed as the composition of two mappings. The first mapping, M µ 1 , is defined to account for the change of radius of the two spheres and it is written in the general separable expression of equation (22) with where x ± 0 = x ± x 0 , R out = 0.45 and, as detailed in section 5.3, x 0 = (1.5, 0) and R ref = 0.116. The radius of the sphere centred at x 0 is defined as R + (µ 1 ) = −0.0372µ 2 1 +0.0968µ 1 + 0.25 so that it takes value 0.116 for µ 1 = −1, 0.25 for µ 1 = 0 and 0.3096 for µ 1 = 1. The radius of the sphere centred at −x 0 is defined in terms of R + (µ 1 ) in such a way that the total volume of the two spheres is maintained, namely (R + ) 3 + (R − ) 3 = 1/32. The piecewise nature of the mapping is illustrated in figure 30, in the vicinity of one of the spheres.
The second mapping, M µ 2 , is defined to account for the change of distance between the spheres and it is written in the general separable expression of equation (22) with where the function d(x) is given by with R int = 0.47 and, as detailed in section 5.3, L = 6.
As illustrated in figure 30 both mappings are defined in a piecewise form. The mappings selected are only C 0 on the artificial interfaces denoted by discontinuous lines in figure 30. Therefore, to facilitate the numerical integration of the terms involving the Jacobian and the adjoint of the mapping, the computational meshes selected are conforming with these interfaces, as it can be observed in the mesh displayed in figure 14.
It is also worth noting that other mappings, with a smooth transition in the artificially created interfaces can be devised. Numerical experiments not reported here for brevity, demonstrate that the piecewise linear mapping described here results in a lower number of integration points required to ensure that errors due to the numerical integration are lower than the interpolation error. However, the choice of a smoother mapping circumvents the need to create meshes conforming with artificially created interfaces. In any case, as stressed in remark 1, this work focuses on the combination of the HDG and PGD formulations and for general geometries the general procedure described in 47 is preferred, rather than the definition of analytical mappings.

D Geometric mapping for the corrugated channel
Similarly to the previous example, the mapping used in the example involving the flow around a sphere in a corrugated channel is designed as the composition of two mappings. The first mapping, M µ 1 , is defined to account for the change of radius of the sphere and it is written in the general separable expression of equation (22) with where R out = 0.4 and R ref = 0.2 and the radius of the sphere, centred at the origin, is defined as R(µ 1 ) = (µ 1 + 2)/10.
The second mapping, M µ 2 , is defined to account for the change of amplitude in the undulatory part of the channel. It only affects the y coordinate and, more precisely, only the definition of f n in equation (57). More precisely, the profile of the channel is given byequation (57) with f n = 1/2 + µ 2 .