Geometrically Parametrised Reduced Order Models for Studying the Hysteresis of the Coanda Effect in Finite Element-based Incompressible Fluid Dynamics

This article presents a general reduced order model (ROM) framework for addressing fluid dynamics problems involving time-dependent geometric parametrisations. The framework integrates Proper Orthogonal Decomposition (POD) and Empirical Cubature Method (ECM) hyper-reduction techniques to effectively approximate incompressible computational fluid dynamics simulations. To demonstrate the applicability of this framework, we investigate the behavior of a planar contraction-expansion channel geometry exhibiting bifurcating solutions known as the Coanda effect. By introducing time-dependent deformations to the channel geometry, we observe hysteresis phenomena in the solution. The paper provides a detailed formulation of the framework, including the stabilised finite elements full order model (FOM) and ROM, with a particular focus on the considerations related to geometric parametrisation. Subsequently, we present the results obtained from the simulations, analysing the solution behavior in a phase space for the fluid velocity at a probe point, considered as the Quantity of Interest (QoI). Through qualitative and quantitative evaluations of the ROMs and hyper-reduced order models (HROMs), we demonstrate their ability to accurately reproduce the complete solution field and the QoI.


Introduction
This study focuses on projection-based reduced order models for fluid dynamics problems with time-dependent geometric parametrisations.We solve the Navier-Stokes equations on a parametrised domain denoted as Ω(µ) ⊂ R d , with a time span [0, T ].Later on, in section 2.1, the specific form of the equations is shown.The geometric parametrisation is defined via a mapping φ(µ) : Ω 0 × P → Ω, such that, given a parameter µ ∈ P ⊂ R p , each point in the original configuration is mapped onto a corresponding point in the deformed one as in Fig. 1.

Bifurcation and the Coanda Effect
The Coanda effect can be described as the tendency of a flow jet to attach to solid surfaces.The term was coined after Henry Coanda, as he patented the first device taking advantage of the effect [59,1].In this paper we adopt as application case a contraction-expansion channel geometry [39,37,8].This case was chosen due to its tendency to exhibit complex dynamics, even at relatively low Reynolds numbers.For one, in this geometry the Coanda effect is present in the form of an asymmetric jet attaching to either the lower or upper walls, depending on the values of the Reynolds number (Re).From the point of view of mathematical nonlinear analysis, this application case presents a supercritical pitchfork bifurcation [3], as shown in Fig. 2. As can be seen, within a certain range of values of Re, a unique solution exists.However, once a critical point is surpassed, it is possible for more than one solution to exist for the same value of Re.For Re ≥ ReSB three solutions coexist.In this paper, our focus lies not in constructing the bifurcation diagram, but rather on studying the application of model reduction to a dynamically morphed contraction-expansion channel operating in a neighbourhood of ReSB The fluid development in a contraction-expansion channel can be described as follows: • For sufficiently small Re, there exists a single perfectly symmetric solution, exhibiting symmetry in both the horizontal and vertical directions.
• As Re increases, the horizontal symmetry of the solution is lost, but the solution remains symmetric about the vertical axis.This results in the formation of a symmetric jet.
• At a critical Reynolds number value, denoted as Re SB (see Fig. 2), the symmetry of the jet is disrupted due to the Coanda Effect.In this scenario, the jet attaches itself to either the upper or lower wall.While the symmetric jet solution is still mathematically possible, it becomes unstable.Numerical simulations or experiments typically yield one of the non-symmetric solutions unless advanced techniques are employed to extract the unstable solution, for instance, when constructing the complete bifurcation diagram [42,41].
• With further increases in Re, the presence of eddies becomes prominent.At a specific critical Reynolds number, denoted as Re H , a Hopf bifurcation occurs, resulting in the emergence of an oscillating solution.However, the discussion of the Hopf bifurcation lies beyond the scope of this work.Interested readers are referred e.g. to [45,51] for further exploration of Hopf bifurcations.
In this investigation, the variation in geometry acts as the sole driving factor for Re.The viscosity and density of the fluid remain constant throughout the study, and once the fluid is initialised, the boundary conditions remain unchanged.In this setting, as the width of the channel narrows, Re increases.

Related Works
The contraction-expansion channel has previously been utilised in research on mitral valve regurgitation disease, as demonstrated in studies such as [43].This disease is characterised by the improper closure of the mitral valve, as depicted in Fig 3 .As a result of this faulty closure, blood is able to flow through a narrow opening, giving rise to the wall-hugging behavior mentioned earlier.Previous numerical studies on the contraction-expansion channel have utilised an affine mapping φ(µ), in conjunction with the spectral element method (SEM) [23].The affine nature of the geometric mapping in that study facilitated a complete online-offline decoupling of the reduced-order models, eliminating the need for additional interpolation and achieving significant speedup factors for the reduced-order models.
In a subsequent study [22], a nonlinear geometric mapping was introduced to account for walls with varying curvature.To effectively evaluate the reduced-order models in this scenario, the discrete empirical interpolation (DEIM) method was employed [7,46].
More recently, the fluid-structure interaction problem has been investigated by incorporating (hyper)elastic walls [30].In that case, the deformation of the walls was computed from their interaction with the flow, rather than being imposed.For the present work, we will adhere to the approach of explicitly imposing the walls deformation.Furthermore, in [41] a non-intrusive neural network framework, known as POD-NN, was employed for efficiently approximating the bifurcating phenomena in a contraction-expansion channel, as well as in other geometries.
In the aforementioned studies, the main objective was to construct the bifurcation diagram.To achieve this goal, steadystate solutions were computed for each parameter variation.This means that even when the geometry was parametrised, the geometric parametrisation remained unchanged over time.The focus was on capturing the steady-state behavior of the system and analyzing the bifurcation phenomena, rather than considering time-dependent variations in the geometric parametrisation.
Time dependent geometry variations for studying wall-hugging phenomena in flow jets have been investigated for example in [2].In that work, both numerical and experimental setups were applied to study the Coanda effect on jets interacting with inclined planes.The parameter to vary was the angle of the plate at the exit of the jet.It was observed that, as the inclination angle of the planes was varied, the jet eventually attached, detached, or re-attached to either one of the walls, depending on whether the inclination angle was being increased or decreased.Thus, resulting on a hysteresis loop.This hysteresis phenomena in fluid problems exhibiting bifurcations had already been reported in classical works [38,31], and more recently in [44,50].

Original Contributions
• Time-dependent geometric parametrisations in a reduced order model context.Our main contributions in this study expand on the existing literature by considering a contraction-expansion channel geometry that un-dergoes time-dependent deformations.By incorporating this dynamically morphed geometry, we are able to observe hysteresis phenomena in the solution, which has been previously documented in works such as [44].
• Independendence of the geometric parametrisation.The framework presented in this work offers broad applicability to a range of geometrically parametrised ROMs.This versatility allows for the framework to be employed in various scenarios.Unlike previous approaches that rely on assumptions about the nature of the mapping φ(µ) to achieve efficient offline-online decoupling, as seen in works like [32,35,46,23], we achieve efficient offlineonline decoupling through the use of the empirical cubature method (ECM) [20,19].Notwithstanding, in this paper we confine our focus on finite element-based ROMs.Consequently, the geometric transformations imposed via the mapping φ(µ) must preserve the integrity of the elements, ensuring the positivity of the Jacobians of the isoparametric transformations.We detail the mappings employed in this study in Sec. 3.For geometric transformations involving more significant deformations incompatible with FEM techniques, alternative ROMs for CutFEM or shifted boundary method are available, as discussed e.g. in [29,28].
• Hyper-reduction technology.We note that the ECM algorithm was designed and has previously been employed for the selection of Gauss points, and the quantity to reduce pertained to internal forces.This paper presents its use in the hyper-reduction of discrete residuals and elements selection.We highlight the similarities of the method with alternative approaches.Furthermore, we emphasises that adopting ECM addresses innately what can potentially be an ill-conditioned problem.

Organisation of this Paper
This paper is structured as follows.In Chapter 2, we present the formulation of the general framework for the fluid problem, including the formulations for the reduced and hyper-reduced order models.Section 3 provides a detailed description of the finite element model utilised in this study, along with the two different geometric mappings employed.One mapping involves an affine transformation with straight walls, while the other employs a nonaffine mapping with curved walls.The obtained results are presented in Section 4. Finally, in Section 5, we discuss the conclusions drawn from this study and provide insights into future research directions.

Formulation
In order to maintain this paper as self-contained as possible, this section presents a concise overview of the stabilised finite element discretisation of the governing equations employed in this study.We then introduce the Galerkin Proper Orthogonal Decomposition (POD-Galerkin) and the ECM hyper-reduction techniques, highlighting the specific considerations related to the geometrically parametrised problem under investigation.We conclude the chapter by providing a summary of the simulation workflow for all the models employed in this study.

Governing Equations
As mentioned earlier in the introduction, we now proceed to explicitly state the fluids problem.The governing equations considered are the standard Differential-Algebraic incompressible Navier-Stokes equations in an Arbitrary Lagrangian Eulerian (ALE) frame of reference where u is the velocity, b are the body forces, σ = −pI + 2ν∇ s u is the Cauchy stress tensor, p is the pressure, ν is the kinematic viscosity, c = u − û is the convective velocity, and û is the so-called mesh velocity which results from the deformation of the domain and is a datum to the fluids problem.The governing equations, taking into account the provided expressions, can be written as follows: The initial conditions g 0 , and Dirichlet and Neumann boundary conditions g D and g N are case-specific and are specified by the corresponding functions, as: with ∂Ω N (µ) ∪ ∂Ω D (µ) = ∂Ω(µ) and ∂Ω(µ) N ∩ ∂Ω(µ) D = ∅.

Full Order Model (FOM)
Let us introduce the functional spaces to pose the weak formulation for the space discretisation as In what follows, we will use the standard L 2 (Ω) inner product notation (•, •) Ω for brevity.The weak form of the governing equations in Eq. 2 is then given by: We define the finite element spaces V h ⊂ V and Q h ⊂ Q.These spaces are spanned respectively by a set of basis functions {ψ i } Nv i=1 and { ψj } Np j=1 , being dimV h = N v , and dimQ h = N p .The finite elements weak form of the governing equations can be expressed as follows: where c h = u h − ûh is the finite elements convective velocity.
To ensure the well-posedness of the problem, the selected spaces for the approximation of pressure and velocity should be compatible and satisfy the Ladyzhenskaya-Babuška-Brezzi (LBB) condition [12].The LBB condition can be stated as follows: inf LBB compatible spaces necessarily comply with however, in this work, we aim to use the same order of interpolation for the finite element approximation of both variables.

Variational Multiscale (VMS) Stabilisation
The VMS approach introduces subgrid spaces V ′ and Q ′ , defined as: where ⊕ denotes the direct sum.The complementary velocity and pressure functions (u ′ , p ′ ) are such that In VMS, the complementary functions are modelled as where the terms τ M and τ C represent respectively the momentum and mass conservation residuals defined as The modelling of the complementary functions with respect to the residuals ensures the consistency of the method, moreover the selection of τ M and τ C depends on the specific type of VMS to apply.For example, considering τ M , τ C = 0 leads to the Galerkin method.In the case of the quasi-static variational multiscale (QSVMS), which we employ in this study, the time evolution of the complementary functions is neglected [10].
To proceed with an algebraic formulation, we consider a spatial tessellation of the domain Ω into N el finite elements, such that e=1 Ω e = Ω; moreover we specify standard finite elements shape functions with compact support and equal degree of interpolation.The velocity and pressure algebraic vectors are obtained by collecting the nodal solutions for both variables as u = (u 1 , u 2 , . . ., u dNn ) T and p = (p 1 , p 2 , . . ., p Nn ) T , where N n is the number of nodes, and d is the number of spatial dimensions.The QSVMS semi-discrete system can then be expressed as follows: where the standard Galerkin finite element terms are given by where A e is the finite elements assembly operator, 1 ≤ i, j ≤ n en , and n en is the number of nodes per element.The remaining terms are related to the VMS stabilisation and are defined as follows:

Time Discretisation
We employ a generalised-α time integration scheme [9].In these schemes, a set of parameters γ, β, α m , α f are defined and used to obtain a linear combination of the terms of the system in Eq. 13 with respect to current and past time steps.Specifically, the Bossak scheme, which is second-order accurate and unconditionally stable, is used in this work.The application of the time integration scheme leads to the fully discrete system: where K eff (d n+ℓ ) ∈ R N ×N is the fully discrete system matrix, N is the total number of degrees of freedom, d n+1 ∈ R N is the finite elements solution vector at time step n + 1, that is The dependence of the matrix K eff (d n+ℓ ) on the solution is due to the nonlinearity in the convective term, the subscript n + ℓ depends on the type of linearisation used.We now say a few words on the linearisation technique employed.

Linearisation of the Discrete System
Given a solution at a time step n, the solution corresponding to time step n + 1 is written in incremental form as where ∆d ∈ R N is a solution increment.We define the residual R : R N × P → R N (highlighting the dependence of all terms of the residual on the geometric parameter µ ∈ P), as together with a fixed-point iteration method of the form where the Jacobian matrix J = ∂R/∂d, and k is the current iterate index.We use a Picard method [12], which amounts to choosing ℓ = 0 for the matrix K eff (d n+ℓ ) in Eq. 19.A Newton-Raphson method would amount to choosing the subindex ℓ = 1, and therefore taking the tangent in the current iterate.

Quantity of Interest (QoI)
It is not uncommon that, rather than the complete solution field d, a quantity of interest is computed through a given operator Γ on the solution field as: Examples of this QoI are the lift and drag coefficients.For the case at hand, such a QoI will be given by the vertical component of the velocity at a selected position indicating the presence of the Coanda effect.

Reduced Order Model (ROM)
Projection-based reduced order models involve constructing an optimal basis that spans a subspace where a highdimensional system can be projected and solved.Various techniques in the literature fall under this classification, and for our research, we have chosen the POD-Galerkin technique.POD, like other projection-based reduced order modeling techniques, comprises two stages: • Offline stage: In this stage, a set of simulations is performed using the computationally expensive FOM, and the resulting solutions are stored in matrix form.These matrices are then analysed to obtain the aforementioned basis.
Ideally, the produced ROMs should not require the evaluation of full dimensional variables.We accomplish this decoupling through a hyper-reduction training.
• Online stage: With the basis and additional hyper-reduction data available, the hyper-reduced order models (HROMs) can be efficiently launched for unexplored parameters at a fraction of the cost associated with the FOMs.
In the subsequent sections, we will delve into the details primarily concerning the offline stage of ROMs and HROMs.

Proper Orthogonal Decomposition (POD-Galerkin)
We define the discrete solution manifold as the set of solution vectors d for all possible values of the parameters vector, that is We employ a linear approximation that, given a reduced order model solution at a time step n, dn ∈ R N , the solution corresponding to the step n + 1 is where ∆q ∈ R N is the reduced solution increment, with N ≤ N (usually N ≪ N is expected) , and Φ = [ϕ 1 , . . ., ϕ N ] ∈ R N ×N , is the reduced basis matrix, obtained by employing the POD method [49,25].
The procedure consists in taking m samples (FOM solutions) of the discrete solution manifold, and store them in a snapshots matrix Here, each of the m samples corresponds to a time step.For this, we consider a function defining a trajectory over the discrete solution manifold with m cases corresponding to the pairs (t, f µ (t)).The specific trajectories used in our investigation are shown later in section 4.2.
Having at one's disposal the snapshots matrix S, we apply the truncated singular value decomposition with a truncation tolerance 0 ≤ ϵ SOL ≤ 1, as S = U N Σ N V T N + E where, The optimal N -dimensional basis [14] is obtained as the truncated matrix of left singular vectors Φ := U N .
Substitution of Eq. 23 into Eq.19, and subsequent projection of the over-determined system of equations onto Φ (Galerkin projection), results in The fixed-point method for solving for the reduced increment ∆q ∈ R N is given as The same operator1 defining the QoI for the FOM, can be used for the reduced order model solution vector as

Hyper-Reduction via Empirical Cubature
Taking into account the finite elements discretisation employed, Eq. 26 can be represented as where L e ∈ {0, 1} e dof ×N is the Boolean operator localising the high dimensional vector of dimension N to the degrees of freedom associated to element e.Consequently, Φ e and de are, respectively, the entries of the basis and reduced solution vector associated to element e, and R e : R e dof s × P → R e dof s is the elemental residual.
Hyper-reduction techniques within finite elements are primarily categorised into two distinct groups, based on their approach to constructing efficient approximations to Eq. 29.On the one hand, there are techniques that focus on approximating (interpolating) the full-dimensional variable R ∈ R N , followed by projecting the approximated vector onto a reduced space.This category, known as the approximate-then-project methods [17], require an orthogonal basis for the nonlinear term, potentially in an unassembled form [57].The process requires employing a greedy algorithm [7] or a QR factorisation with pivoting [13] to identify critical sampling points, often referred to as "magic points" [53].
On the other hand, there are strategies that attempt to approximate directly the projected vector Φ T R ∈ R N .These techniques are classified as project-then-approximate [17] methods.Here, the approximation of Eq. 29 is performed not by sampling the full set of elements, but rather a selected subset E ⊂ {1, 2, . . ., N el }.Each considered projected elemental contribution is in turn adjusted using a corresponding weight factor ω e .
We discussed this classification in detail in [20], where we used the concepts of nodal vector approximation approaches and integral approximation approaches.The methodology adopted in the current paper can be classified as an integral approximation or a project-then-approximate method.Our preference for this type of methods is motivated by their demonstrated stability and robustness, which are discussed in more depth in the subsequent paragraphs.
Having summarised the classification of the hyper-reduction techniques, we proceed to detail the necessary steps to obtain the reduced set of elements and weights.We begin by storing the elemental contributions in Eq. 29 for each of the m studied cases.Let the projected residual for element e for parameter case i (here µ i = f µ (t i ), see Eq. 24) be defined as We construct then the matrix of projected residuals for all N el elements and m studied parameters, as The exact assembly of the projected residuals, for the m studied parameters variations, is given as where The optimisation problem for finding the reduced set of elements E and corresponding weights ω is then posed as where ζ ∈ R N el is a sparse vector with non-zero values ω at indices E, ∥•∥ 0 represents the zero pseudo norm (counting the number of non-zero entries of its argument), ϵ is a user-defined tolerance, and we use the symbol • ⪰ • to represent inequality with respect to the non-negative orthant (this constraint requires the non-negative entries of ζ to be positive).This non-negative constraint was considered critical for maintaining the positive definiteness of the involved operators in works such as [15,20].Moreover, a proof that solutions to the optimisation problem in Eq. 33 are stable was provided in [16], showing that such solutions preserve the Lagrangian structure of the problem at HROM level.
Given the optimisation problem, the question arises as to which algorithm should be used to obtain a solution.It is widely recognized [6] that the problem, as formulated in Eq. 33, is NP-hard, hence computationally intractable.Therefore, one must resort to suboptimal greedy heuristics or convex relaxation.Among the alternatives is the Energy-Conserving Sampling and Weighting (ECSW) algorithm [15,16], which employs a non-negative least squares (NNLS) procedure with an early stopping criterion to account for the user-defined tolerance.Another alternative involves modifying the cost function to include the ∥•∥ 1 norm, thus convexifying the problem and making it amenable to established convex optimisation algorithms, as suggested e.g. in [40,60].
The ECM algorithm, that we employ in this study, addresses a modified version of the optimisation problem by initially applying a truncated SVD, with truncation tolerance ϵ RES to matrix S r , as where We use the truncated matrix of right singular vectors G to define a vector containing the assembly of the modes of the projected residuals (in contrast to the exact assembly of the projected residuals in Eq. 32), as In this way, the optimisation problem tackled by the ECM algorithm is By operating on a basis for the row space of S r , the ECM algorithm circumvents the need to employ a NNLS algorithm and instead uses a standard least-squares algorithm in every greedy iteration.Admittedly, the computational burden is now passed to the SVD algorithm, which could be required to tackle large matrices.Fortunately, efficient SVD algorithms have been proposed; for example, a sequential-randomized SVD for dealing with matrices larger than the available memory on a local PC was presented in [26].Parallel algorithms are also available for the SVD, for example leveraging the talland-skinny matrices involved [58].
An additional advantage of employing the ECM is that it readily incorporates a treatment for ill-posed problems, consisting of assembled projected residuals close to or numerically zero, i.e., when S r 1 ≈ 0 [54].This condition was discussed in our paper [20] for self-equilibrium problems, where the equivalent quantity for the assembled projected residuals is exactly zero by construction.As a result of this prevision, the application of the ECM algorithm leads to the selection of β + 1 elements, one more than the number of rows in G, as Such solution exactly integrates the assembled modes of the projected residuals as A discussion on the relation between the truncation tolerance ϵ RES in Eq. 34, and the user-defined tolerance ϵ in the original optimisation problem in Eq. 33 has been recently presented in [26].One can readily check that a solution obtained by the ECM algorithm also complies with Indeed, when using the ECM algorithm, the truncation tolerance for the SVD of the projected residuals, ϵ RES , is the parameter dictating the sparsity of the reduced mesh, and in turn the precision of the HROM integration with respect to the ROM.If no truncation is performed, the number of selected elements will equal Rank(S r )+1.
Having at one's diposal the HROM elements and weights, the reduced terms in Eq. 27 are assembled by looping over the reduced set of elements multiplied by the corresponding elemental weight, as

Global Workflow
Having introduced all the components, we now straightforwardly assemble them in this section as a summary.The following box presents the summarised workflow for both the offline and online stages of the POD-Galerkin with ECM Hyper-reduction.

Offline stage 1. Define the m training scenarios {µ
Launch the corresponding FOM simulations and build the snapshots matrix S ∈ R N ×m 3. Compute the SVD of S, with truncation tolerance ϵ SOL to get the basis Φ ∈ R N ×N 4. Construct the snapshots of residuals projected S r ∈ R N •m×N el 5. Compute the SVD of S r , with truncation tolerance ϵ RES , to get the basis for its rowspace G ∈ R β×N el 6. Obtain the set of elements and positive weights using the ECM algorithm [19] (E, ω) ← ECM(G, 1) 2. Construct and solve the reduced system by looping over the subset of elements E multiplied by their corresponding weight ω

Box 2.1: Workflow of POD-Galerkin with ECM hyper-reduction
The m training scenarios described in step 1 of the offline stage in Box 2.1 are defined in this work using a specific function f µ defining a training trajectory, as formulated in Eq. 24.It is possible to employ more than one training trajectory; furthermore, incorporating other parameter variations can be achieved without significant complications.Similarly, in this work, the parameters introduced in step 1 of the online stage, denote testing trajectories.These trajectories are distinct from those used during the training stage.
Alg. 1 shows the speudocode for running the µ−geometrically parametrised simulation, where the steps to be followed by either a FOM, ROM or HROM simulation are delineated.The main difference among these models occurs in point 7.
Algorithm 1: Pseudo code for the geometrically parametrised simulations 1 Input: final time T , time increment size ∆t, trajectory function f µ (t) : t → µ, for a ROM also the basis Φ ∈ R N ×N , and for an HROM, both the basis Φ ∈ R N ×N and the HROM elements and weights E, ω 2 Output: State variable u n ∈ R N for all time steps (S) and/or QoI z n ∈ R α for all time steps (Z) update the geometric parameter: compute the mesh velocity ûh ← g mesh (Ω n+1 ) get the solution increment ∆d n+1 by • FOM: solve until convergence system in Eq. 20 • ROM: solve until convergence system in Eq. 27 • HROM: solve until convergence system in Eq. 27, construct the system as in Eqs.40 and 41 8: It is worth noting that for this work, step 4 of Box 2.1 was performed by running a ROM following Alg. 1, while storing the projected residuals as outlined in Section 2.3.2.This way of performing the HROM training can be expensive (in particular in terms of memory).In other works, the authors have chosen to project the readily available snapshots matrices onto the column space of the basis as d projected = ΦΦ T d, to finally obtain the residuals with respect to these projected snapshots.Our decision to run Alg. 1 to obtain the residuals projected is completely consistent with the theory presented in section 2.3.2.

Model Description
The plane contraction-expansion channel used in this study is shown in Fig. 4.This geometry serves as the base or reference configuration Ω 0 .To effectively morph this base geometry into a deformed configuration Ω, as depicted in Fig. 1, we employ two distinct geometric mappings.This approach is employed to showcase the framework's versatility, irrespective of the nature of the mappings.Each geometric mapping depends on a single scalar geometric parameter µ ∈ R, and are in turn an affine mapping denoted as φ AFFINE and a nonaffine mapping that combines Free Form Deformation (FFD) and Radial Basis Functions (RBF) denoted as φ FFD + RBF .
The affine mapping φ AFFINE changes the narrowing width w c in a symmetric way while preserving the walls straight.On the other hand, the nonaffine mapping φ FFD + RBF allows to obtain deformed configurations that undergo more complex transformations and lead to curved walls of the narrowing.

Affine Mapping
For the affine mapping, the only parameter that dictates the shape of the geometry is the narrowing width µ = w c ∈ [0.1, 2.9].Following [23], the affine mapping is defined by decomposing the original domain Ω 0 into N dom non-overlapping subdomains as, and defining for each subdomain For the geometry at hand, five distinct regions can be identified, as shown in Fig. 5.For each of them, a particular transformation matrix allows to map the point in the original configuration to the deformed configuration, using as parameter the narrowing width w c .Eq. 45 shows the transformation matrices corresponding to each of the coloured regions.

Nonaffine Mapping
We perform the nonaffine mapping in two stages.In the first stage, an FFD technique is applied to the boundary points of the 2D geometry.The second stage involves moving the points in the interior of the geometry.This is accomplished using an RBF interpolator.
In FFD, the geometry to morph is surrounded by a box of control points, as can be seen in Fig. 6.Let P 0 represent the set containing the original position of the control points defining the bounding box.For our particular case, the set of deformed control points' positions P (µ) depends on a single scalar µ ∈ [0, 1].The deformed position x ∈ ∂Ω of a point x 0 ∈ ∂Ω 0 is then given as, For information on the specific expression of the mappings and blending functions used, the interested reader is directed to [48,34,56].We employ the implementation of FFD from the Python library PyGeM [55].After the deformed boundary points have been computed, the deformed position of the points in the interior of the geometry x ∈ Ω can be obtained by applying an RBF interpolator.
Let xi represent the i-th undeformed boundary mesh point (whose deformed position is obtained via the FFD mapping).The deformed position x ∈ Ω of a point x 0 ∈ Ω 0 is then given as, where N b is the number of points in the deformed boundary, ϕ(•) is a given basis function depending on the distance from the desired point in the interior to a mesh point in the boundary, and q(x 0 ) is a polynomial.The coefficients β i and the polynomial q(x 0 ), are obtained as a function of the deformed boundary points by imposing the interpolation conditions Further details about RBF interpolation can be found in [53,5].Successive application of the FFD mapping for all mesh points on the boundary of the geometry, followed by RBF parameters computation and application for all mesh points in the interior, completely defines the mapping (49)

Mesh
The two geometric mappings require different control surfaces with a corresponding label to exist.For example, the affine mapping requires to clearly identify the coloured regions in Fig. 5 to apply the respective transformations on them.Therefore, a conforming mesh fitted to each of them was generated.On the other hand, the nonaffine mapping only required to identify the upper and lower rectangles defining the narrowing.Moreover, while the nature of the affine mapping allows to preserve the positive Jacobians of the isoparametric finite elements under relatively large deformations, the nonaffine mapping can lead to distorted elements which intersect, or completely swap, whenever deformations are too large.These differences on the mappings forced us to use a different mesh for each of them, while keeping the same element size in both.We have used unstructured meshes in both cases because they are known to work better for the specific example under consideration [45].The information of the meshes used2 can be seen in Table 1, while Fig. 7 show the generated meshes.In both figures, the presence of the probe point p * allows to monitor the evolution of the QoI, and therefore the occurrence of the Coanda effect.As the geometry morphs in time following a circular trajectory, a hysteresis plot can be created.The position of the probe point is selected to be located at a node in the corresponding finite element mesh that is close to the y-direction centre line, and slightly behind the narrowing.The generation of the meshes was performed using the pre-and post-processor GiD [36].The general purpose software KratosMultiphysics [18] was used for launching the FOM, ROM and HROM simulations, by taking advantage of the KratosRomApplication.In Kratos, the finite elements discretisation is closely linked to the physics to be simulated.For the physics in these simulations, Navier-Stokes in an Arbitrary Lagrangian Eulerian framework, the only available option was to use P1-P1 triangular elements.

Boundary and Initial Conditions
Both models prescribe the initial and boundary conditions shown next As can be seen, the maximum horizontal velocity at the inlet is set to 3.Moreover, the inlet velocity has been allowed to undergo an initialisation period for one second following a smooth path.Finally, both upper and lower walls prescribe non-slip conditions, and we are not considering body forces.

Material Properties
The values of the viscosity and density for the fluid have been selected so that, given the distortion allowed by the finite element discretisations under the respective geometric mappings, the results were qualitatively comparable to the ones reported in similar works, e.g.[45,39,43,22,23,21].Moreover, the material properties are held constant for all simulations.
For both models, we used the values of dynamic viscosity and density reported in the following table .fluid property value dynamic viscosity μ 0.1 density ρ 1.0

Results
In this section, we begin by providing a qualitative overview of the FOM.We present the different geometrical configurations and showcase the observed Coanda effect along with their corresponding hysteresis loops.Additionally, we offer insights into the training process of the ROMs and HROMs.We discuss the different trajectories and present the results obtained from all models providing both a qualitative and a quantitative comparison including error and speedup factors for each of the two examples studied, which are indented to highlight the performance of ROMs and HROMs in challenging scenarios.

FOM Results
Before delving into the study of the hysteresis of the Coanda effect for the contraction-expansion channel under scrutiny, we first examine the behavior of the FOM and how variations in the geometry impact the fluid solution.
Fig. 8 illustrates the flow field for extreme values of the geometric parameter µ for both the affine geometric mapping (ranging from 0.1 to 2.9) and the nonaffine geometric mapping (ranging from 0 to 1) in the contraction-expansion channel.
For the case where the narrowing width w c is set to 1 (Figures 8b and 8d), both mappings yield qualitatively similar solutions, characterized by a horizontally symmetric jet.However, as the narrowing width decreases, the Coanda effect becomes present, and small eddies appear on the asymmetric jets (Figures 8c and 8f).Fig. 8c exhibits a complex flow pattern.Notice the relatively small value of the narrowing width made possible by the affine geometric mapping thanks to its bijectivity.On the other hand, in Fig. 8f, the maximum distortion allowed by employing the nonaffine geometric mapping produces a narrowing width of w c = 0.422.Although this value is larger than   In Fig. 9, we show the possible solutions obtained for the same values of the geometric parameter µ for both mappings.Advanced methods can be used for obtaining one or the other of the stable solutions, e.g.[42].For our specific case, we have observed that given a specific trajectory (earlier mentioned to be defined by f µ ), that is, given an initial value of the parameter µ, and then deforming the geometry in a closed loop, one of the stable solutions is obtained consistently, and a hysteresis loop can be observed.

Trajectories
Here we highlight the different values of the parameter µ imposed to both models, through the description of the trajectory functions f µ : t → µ.
In this work, we are working with three distinct trajectories.The first two trajectories are straightforward and represent a contraction followed by an expansion, and an expansion followed by a contraction respectively.Both models are therefore subjected to a morphing going from an initial geometry Ω 0 , going to a maximally deformed geometry and from there back to the initial geometric configuration.These first two trajectories are therefore Trajectory 1: Contract then expand Trajectory 2: Expand then contract The third trajectory is a more complicated function defined as follows

Trajectory 3: Trigonometric function
For all of the trajectories, we keep the initial geometry static for a period of time T 1 to allow the flow to develop.Moreover, for trajectories 1 and 2, we treat differently time instances T 2 and T 3 .One the one hand, for the nonaffine mapping we keep the value of the parameter µ static during enough time for the Coanda effect to fully develop.On the other hand, in the case of the affine mapping the asymmetric jet is already strongly present in the solution by the time the maximum deformation is reached, therefore, for this mapping we immediately start the progressive deformation back to the initial configuration.
A graphical representation of these functions applied to each of the two models is shown in Figs. 10, 11, and 12.
where the snapshots matrix S ∈ R N ×m , and the basis matrix Φ ∈ R N ×N , as described in Sec.2.3.1.
Moreover, for each of the two models and for each of the four ROMs, we construct HROMs employing four different truncation tolerances for the SVD of projected residuals, namely ϵ RES = {1e − 3, 1e − 4, 1e − 5, 1e − 6}, such that where the snapshots matrix of projected residuals S r ∈ R N •m×N el , and the basis for its rowspace G ∈ R β×N el .Finally, as described in 2.3.2, the ECM algorithm is used to obtain the HROM elements and weights, as Therefore, for each of the two examples, 8 ROMs and 32 HROMs have been constructed, half of them for the model employing the affine geometric mapping and the other half corresponding to the model employing the nonaffine mapping.
To asses these ROMs and HROMs, we employ three strategies: • Qualitative assessment through phase space plot Throughout this section we use a phase space plot for the QoI, which as mentioned before is the horizontal velocity at the probe point p * (as shown in Fig. 7).The plot consists then on a comparison between the horizontal velocity v * y , against the narrowing width w c .This plot straightforwardly provides insight on the performance of the ROMs for reconstructing the behavior of the FOMs, as well as of the performance of the HROMs for reconstructing the respective ROMs and FOMs.

• Quantitative assessment through measure of errors
To quantitatively evaluate the results obtained for the different models, we define a percentage error operator as follows: Here, a and b can be either matrices or column vectors.In our case, we will use this operator to compare the snapshots matrices of the solution fields for the FOM, ROM, and HROM, denoted respectively as S FOM , S ROM , and S HROM .As well as for the QoI, denoted as v * y FOM , v * y ROM , and v * y HROM for each of the considered models.

• Speedup
In order to quantify the speedup factor, a comparison was conducted between the ROM and the HROM in relation to the FOM.This analysis considered the time required for the construction and solution of the linear system of equations.The time necessary for visualising the complete solution field was not taken into account.
To evaluate the speedup factor, we introduced a speedup operator denoted as s, as follows: where T a represents the time required for either the ROM or the HROM, while T FOM corresponds to the time taken by the FOM.The output of this operator indicates the number of times the reduced order model is faster in comparison to the FOM.

Example 1
For the first example, Trajectory 1 serves as the training trajectory, while Trajectory 2 is taken as the testing trajectory.Fig. 13 displays the QoI phase space plot for both models and for both trajectories.The expectation is that ROMs whose modes favor one branch in the bifurcation encounter difficulties accurately capturing the behavior of the opposite branch.Consequently, this testing approach offers valuable insights into the performance and limitations of the ROMs and HROMs in capturing the full range of behavior of the solution.14 illustrates the singular values of the solution snapshots and the corresponding number of POD modes required by both models.These results are presented for various truncation tolerances ϵ SOL .Notably, the nonaffine geometric mapping consistently demands a greater number of modes to achieve the same tolerance compared to the affine mapping.We generated four ROMs for each geometric mapping, accounting for the respective number of POD modes.We now check the performance of the ROMs for reconstructing the training trajectory, both qualitatively and quantitatively.Fig. 15 shows the QoI phase space plot of the FOMs and ROMs, demonstrating the consistent behaviour of the ROMs: as the truncation tolerance ϵ SOL decreases, the discrepancy between the QoI of the ROMs and that of the FOMs also decreases.Additionally, Fig. 16 illustrates the percentage errors for the reconstruction of the QoI and the complete solution fields, quantifying the observed trend in the phase space plots.We now turn our focus to the performance of the ROMs for reconstructing the testing trajectory.Fig. 17 presents the QoI phase space plots for FOMs and ROMs.We note that none of the ROMs for the affine mapping accurately selects the "correct" branch of the bifurcation, while the ROMs for the nonaffine mapping progressively converge towards the FOM solution.18 shows the percentage errors incurred by both models during the reconstruction of the QoI and the complete solution fields of the respective FOM.We note that the errors in the QoI are significantly larger than those in the complete solution fields.For the affine mapping, the errors in the QoI of the ROMs reach around 145%, whereas the errors in the complete solution fields are approximately 9%.In the case of the nonaffine mapping, the ROMs result in smaller errors in the QoI, except for the model with a truncation tolerance ϵ SOL = 1e − 3, which exhibits the opposite branch in the bifurcation compared to the FOM.

HROM
Since the ECM algorithm computes rules featuring one element more than the number of components in the basis for the rowspace of the matrix of projected residuals (See Sec.2.3.2), the decay profile of its singular values provides an indication of the potential to hyper-reduce the case at hand.Fig. 19 shows the singular values decay profile of S r and the number of selected elements for of the HROMs created for the affine geometric mapping.Moreover, Fig. 20 shows it for the nonaffine mapping.In both cases, we can observe that the larger the value of the truncation tolerance ϵ SOL , the more pronounced the decay profile.This affects the required number of components of the SVD of S r to comply with the four truncation tolerances ϵ RES , and therefore directly dictates the number of selected elements by the ECM algorithm.For each geometric mapping, we have constructed 16 HROMs containing the number of elements specified on letter (b) in the above figures.Appendix 6 shows the location of the selected elements in the meshes employed.We now proceed to check their performance to reconstruct the training trajectories.Fig. 21 shows the QoI phase space plot for each of the ROMs against their corresponding HROMs for the affine mapping.We recall that the error committed by the HROM is bounded by the corresponding ROM used for its construction, however, in these plots we have also included the FOM with a ghosted line for reference.We immediately observe that some of the models are numerically unstable.In the plots, exploding solutions were also ghosted to allow visualisation of the rest.In particular, all of the HROMs constructed with a truncation tolerance of ϵ RES = 1e − 3, as well as the HROM constructed with the combination ϵ SOL = 1e − 3, ϵ RES = 1e − 4 were not numerically stable.The rest of the HROMs behave as expected, approaching the corresponding ROM solution, as the tolerance ϵ RES gets smaller.Turning our attention to the HROMs associated with the nonaffine geometric mapping, we evaluate their performance for reconstructing the training trajectory by first observing Fig. 23.This figure shows the QoI phase space plots of the ROMs against their respective HROMs.We observe that all of the models are numerically stable, but some of the HROMs choose the opposite stable branch in the bifurcation with respect to the ROM.Having assessed the performance of the HROMs for reproducing the training trajectory, we now turn our attention to their capability to handle the testing trajectory.
We begin by examining the HROMs employing the affine geometric mapping.Fig. 25 presents the QoI phase space plot of each of the four ROMs against their corresponding HROMs.As in the case of the training trajectory, all of the HROMs constructed with a tolerance ϵ RES = 1e − 3 are numerically unstable, and as ϵ RES gets smaller, the HROMs approach their corresponding ROMs.For this case, however, the FOM selects the opposite stable branch in the bifurcation diagram with respect to both ROMs and HROMs.Fig. 26 provides quantification of these trends, for example although the errors of the HROMs with respect to their ROMs can be less than 1% in the QoI, the errors with respect to the FOM are no less than 145%.Turning our attention to the performance of the HROMs for reconstructing the testing trajectory with the nonaffine geometric mapping, Fig. 27 displays the QoI phase space plot of the ROMs against their corresponding HROMs.Remarkably, all HROMs are numerically stable, progressively approaching their corresponding ROMs as the tolerance ϵ RES decreases.Fig. 28 offers quantification of the errors for QoI and solution fields of the HROMs against ROMs and FOM.Fig. 29 shows the speedup factors obtained for all models with respect to the FOMs of each geometric mapping, as computed using Eq.56.Notably, the ROMs without hyper-reduction achieve a speedup of less than an order of magnitude.In contrast, the HROMs exhibit significantly larger speedup factors.It is important to exercise caution, considering that the HROMs constructed considering the less stringent truncation tolerances resulted invariably in high errors.The results presented for this first example highlight the challenges of accurately capturing the Coanda effect and its hysteresis behavior using reduced order models whose testing trajectories require capturing the behavior of the opposite stable branch in the bifurcation, with respect of the information contained in the training trajectory.While in the case of the affine geometric mapping, none of four ROMs was capable of accurately reproducing the FOM trend in the Coanda effect for the testing trajectory, ROMs for the nonaffine geometric mapping were able to capture the main patterns of this unseeen brach when sufficiently decreasing the truncation tolerances of the solution and residuals.In this case, although none of the ROMs perfectly matched the FOM solution, they progressively approached it.
To delve deeper on the potential of ROMs and HROMs to deal with the characterisation of the Coanda effect in this geometry, we proceed to consider both Trajectory 1 and Trajectory 2 as training trajectories in Example 2.

Example 2
In the second example, we consider the use of two training trajectories: Trajectory 1 and Trajectory 2. The presence of jets attaching to both, the upper and lower walls in these training trajectories anticipates that the ROMs and HROMs will be robust.We subject all models to Trajectory 3 as the testing trajectory and asses their performance to check this hypothesis.

ROM
Fig. 30 illustrates the singular values of the solution snapshots and the corresponding number of POD modes required by both models.As in the previous example, the nonaffine geometric mapping requires a greater number of modes to achieve the same tolerance compared to the affine mapping.We generated four ROMs for each geometric mapping, accounting for the respective number of POD modes.We asses the capability of the constructed ROMs to replicate the training trajectories, starting with Trajectory 2. Fig. 31 and Fig. 32 show the QoI phase space plots and the percentage errors for the ROMs against the FOMs for both geometric mappings.The ROMs for both geometric mappings exhibit the expected behaviour that, as the truncation tolerance ϵ SOL decreases, the discrepancy between the QoI of the ROMs and that of the FOMs also decreases.36 show the singular values decay profile of the matrices of projected residuals, and the count of selected elements by the ECM algorithm across the 32 constructed HROMs.Appendix 6 shows the location of the selected elements in the meshes employed.As anticipated, the number of HROMs surpasses those in the preceding example.This is a consequence of the richer POD bases, and the subsequent slower decay of singular values of matrices S r (ϵ SOL ).In this example, the results obtained for both the ROMs and HROMs in reproducing the training trajectories are consistent that is, as the tolerances are reduced, they converge towards their respective ROMs and FOM.Moreover, in the case of Trajectory 1, the results are comparable to those from the previous example for the same trajectory.This similarity is anticipated, given that the ROMs and HROMs are subjected in both cases to the trajectories used for their construction.To thoroughly assess their capabilities, we proceed to subject the ROMs and HROMs described in this example to Trajectory 3.

Trajectory 3
We now examine the performance of the ROMs for each geometric mapping in reproducing Trajectory 3, considered as the testing trajectory for this example.The results are depicted in Fig. 45 and Fig. 46.
For the affine geometric mapping, Trajectory 3 induces solutions attaching to both the upper and lower walls, posing a challenge for the ROMs to capture accurately.Notably, the ROM with the smallest error in the QoI is not associated with the smallest truncation tolerance ϵ SOL , but rather to the second smallest.However, the percentage error in the solution field is the same for both, at 10.1%.
In contrast, the ROMs for the nonaffine geometric mapping effectively capture the behavior of the solutions, demonstrating errors in both the QoI and the solution field of less than 1% for the smallest truncation tolerances.In this example, we examined the performance of ROMs and HROMs created using a training set that includes data from both stable branches of the bifurcation.Specifically, with snapshots containing a jet attaching to the upper and lower walls.The results showed that these models can accurately represent not only the two training trajectories but also a more complex testing trajectory.

Summary of Results
As a way of summary, we can highlight three main points: • On the richness of the training sets.In order to develop robust ROMs and HROMs capable of adapting to unseen trajectories, it is essential to include snapshots from both stable branches of the bifurcation in the training sets.This necessity becomes evident when comparing the results from both examples.In example 1, the testing trajectories resulted in significant accuracy issues, because the FOM for the testing trajectory presented solutions on the opposite branch with respect to the training trajectory.In contrast, Example 2, which included both stable branches in its training trajectories, demonstrated improved accuracy for unseen testing trajectories.
• On the truncation tolerance of the solutions ϵ SOL .The truncation tolerance of the snapshots matrix of solutions sets an upper bound on the achievable accuracy of the ROM, and subsequently, its derived HROMs.Opting for larger truncation tolerances does not inherently result in numerically unstable models, as observed in all studied ROMs, where numerical stability was maintained regardless of the incurred errors.However, in scenarios involving both stable branches of the bifurcation, and therefore not limited by this fact to preferring a given solution, such as in Example 2, selecting excessively large values for ϵ SOL caused the ROMs to choose the opposite stable branch compared to the FOM.While this deviation is still physically permissible, it may be significant enough to render some models unsuitable depending on the specific application.
• On the truncation tolerance of the residuals ϵ RES .HROMs constructed with a truncation tolerance of ϵ RES = 1e − 3 exhibited considerable deviations from their corresponding ROMs, to the extent of becoming numerically unstable.This observation holds true, especially in the context of Example 1, but it remains applicable even in scenarios where the POD bases were enriched with information from both stable branches of the bifurcation, as in Example 2. Opting for a more restrictive tolerance, such as ϵ RES ≤ 1e − 4, appears necessary to ensure that the HROMs provide accurate solutions with respect to their ROMs.
Based on these observations, the adoption of the presented HROM framework for characterising the hysteresis of the Coanda effect appears to involve a trade-off between accuracy and efficiency.Fast HROMs, while capable of capturing the overall solution trend, may struggle to select the same branch as the FOM, even when such a solution is achievable the POD basis.On the other hand, more accurate models that reproduce intricate solutions can be obtained by choosing sufficiently stringent truncation tolerances, albeit resulting in smaller speedup factors.
Finally, we acknowledge that we detected a significant increase in assembly time in our implementation as the number of POD modes increased.In scenarios where a relatively large number of modes are involved, we deviated from the element-by-element approach outlined in Section 2.3 for assembling the system of equations.Instead, for these cases we adopted a "global" approach.This entailed assembling the sparse system matrix in a manner similar to the FOM, followed by sparse-dense, and posterior dense-dense matrix product with the complete basis matrix Φ.We observed that the element-by-element formulations consistently outperformed this global approach for cases considering less than 60 POD modes, and a reduced mesh containing less than one fifth of the total elements.The speedup factors reported at the end of each example show the superior formulation between the two options.

Conclusions and Perspectives
In this paper, our focus was on investigating a general ROM framework for addressing fluid dynamics problems with time-dependent geometric parametrisations.This framework encompasses the utilisation of two powerful techniques: POD-Galerkin together with ECM hyperreduction.By employing these techniques, we aimed to effectively capture the intricate fluid behavior inherent in the contraction-expansion channel geometry.While this geometry offers a relatively straightforward setting, it still presents complex fluid dynamics phenomena, such as a bifurcating solution known as Coanda effect.
By utilising ROMs and HROMs, we have successfully constructed accurate models capable of capturing both the training trajectories, which represent a specific deformation of the geometry over time, and challenging testing trajectories, which either trigger the opposite branch of the bifurcation compared to the training trajectory, or introduce a more complicated deformation sequence.
We have analysed the solution behavior in a phase space, specifically focusing on a QoI, which is the velocity in the y-direction at a probe point.This QoI allows for the detection and characterisation of phenomena such as the Coanda effect and its hysteresis.By qualitatively assessing the outputs of the ROMs and HROMs in this phase space plot, we gain insights into the performance of the models.Additionally, quantitative evaluations have been conducted to assess the accuracy of the complete solution field and the QoI.
As discussed in Sec.4.6, the behaviour of the ROMs and HROMs suggest a trade-off between accuracy and computational speedup.The HROM models exhibit significant speedups, while still providing physically acceptable and bounded solutions.However, these models incur relatively large errors in reproducing the complete solution field and the QoI of FOMs, particularly for testing trajectories.Despite these errors, these models can still be useful in applications where the detection of the Coanda effect is crucial, even if the selected bifurcation branch is incorrect.This is evident, for instance, in Example 2, where the HROMs successfully capture the overall trend of the solution for unseen testing trajectories, despite the fact that the obtained solution does not match the FOM solution.For more accurate results, HROMs offering moderate speedups while maintaining low errors can be employed.For applications demanding very high levels of accuracy for unseen trajectories, along with very high speedup factors, the incorporation of nonlinear ROM techniques could be considered, as outlined in Future Work.

Future Work
There are several promising avenues for further advancing this research.
Firstly, it is important to acknowledge that our study focused on a single parameter variation within a mesh comprising a relatively small number of elements, which was suitable for our academic objectives.However, in scenarios where higher resolution and increased accuracy are required, it becomes imperative to employ a larger number of elements in the base model.Additionally, an effective reduced order model should be trained by exploring a multidimensional parameter space.Addressing these considerations necessitates addressing certain challenges within the framework presented in this paper.Specifically, the launch and analysis of simulations, as well as the management of the generated matrices via singular value decomposition, become computationally demanding on a single machine.
As mentioned in Section 2.3.2, the size of the snapshots matrices increases with the number of elements and the number of POD modes.To alleviate this challenge, we are actively engaged in the development of parallelisation techniques for the entire workflow.This includes parallel simulation orchestration, efficient data management strategies, and the implementation of parallel algorithms for computing the singular value decomposition.These parallelisation efforts aim to significantly enhance the computational efficiency and scalability of the training process, enabling the exploration of larger parameter spaces and higher fidelity models.Furthermore, simulations involving qualitatively distinct solutions, such as the ones demonstrated in this paper, often require a large number of modes to capture the intricate behavior of the solutions.To address this challenge, we are exploring alternative strategies.One approach involves utilising multiple piece-wise linear bases that effectively capture the specific behavior in the vicinity of a particular region in the parametric space, as demonstrated in [21,24].Additionally, we are investigating the utilisation of nonlinear manifolds, particularly quadratic approximations as proposed in [27,4], to mitigate the requirement of a high number of modes.Lastly, we are actively exploring the application of autoencoder neural networks as a form of generic manifold Galerkin approximation [33,47].We anticipate that the results of these advancements will be reported in subsequent papers, expanding upon the findings presented in this study.

Figure 1 :
Figure 1: Transformation of a domain Ω 0 into Ω via a mapping φ

Figure 2 :
Figure 2: Supercritical pitchfork bifurcation diagram for a contraction-expansion channel.For Re < ReSB there exists a unique solution.For Re ≥ ReSB three solutions coexist.In this paper, our focus lies not in constructing the bifurcation diagram, but rather on studying the application of model reduction to a dynamically morphed contraction-expansion channel operating in a neighbourhood of ReSB

Figure 3 :
Figure 3: Mitral valve of the human heart and the contraction-expansion channel model used in this investigation.Contraction-expansion channel models can help in the study of mitral valve regurgitation disease

==
∇ s ψ j , ∇ s ψ i Ω e the diffusion matrix , C(c) = A e C e C e ij = − ψ j , c h • ∇ψ i Ω e the convection matrix , − ψj , ∇ • ψ i Ω e the pressure/divergence matrix , F = A e F e F e i = (b, ψ i ) Ω e the forcing term ,

Figure 4 :
Figure 4: Base model of contraction-expansion channel

Figure 5 :
Figure 5: Reference subdomains for affine mapping

Figure 6 :
Figure 6: FFD part of the nonaffine mapping.For the second part of the nonaffine mapping, the interior points are moved according to an RBF interpolator (a) Mesh for mapping φAFFINE (b) Mesh for mapping φFFD+RBF

Figure 8 :
Figure 8: FOM solutions for different deformed configurations

Figure 13 :
Figure 13: QoI phase space plots depicting the impact of the two trajectories in the FOMs for both geometric mappings.a) When subjected to Trajectory 1, the FOM for the affine mapping produces a solution with a jet attaching to the lower wall of the channel, while applying Trajectory 1 to the FOM of the nonaffine mapping results in a jet attaching to the upper wall.b) When subjected to Trajectory 2, the FOM for the affine geometric mapping produces a jet attaching, first to the upper wall and subsequently to the lower wall.On the other hand, when applying Trajectory 2 to the FOM for the nonaffine mapping, the jet produced attaches only to the lower wall, before going back to the symmetric position

Figure 14 :
Figure 14: Singular values of solution snapshots and the corresponding number of POD modes for both geometric mappings under Trajectory 1

Figure 15 :Figure 16 :
Figure 15: QoI phase space plot for the FOM against various ROMs for the training trajectory (Trajectory 1) for both geometric mappings

Figure 17 :
Figure 17: QoI phase space plot for the FOM against various ROMs for the testing trajectory (Trajectory 2) for both geometric mappings

Fig.
Fig.18shows the percentage errors incurred by both models during the reconstruction of the QoI and the complete solution fields of the respective FOM.We note that the errors in the QoI are significantly larger than those in the complete solution fields.For the affine mapping, the errors in the QoI of the ROMs reach around 145%, whereas the errors in the complete solution fields are approximately 9%.In the case of the nonaffine mapping, the ROMs result in smaller errors in the QoI, except for the model with a truncation tolerance ϵ SOL = 1e − 3, which exhibits the opposite branch in the bifurcation compared to the FOM.

Figure 18 :
Figure 18: Percentage error on QoI and solution field of FOM against various ROMs for the testing trajectory (Trajectory 2) for both geometric mappings

Figure 19 :
Figure 19: Singular values decay profile of matrices of projected residuals, and number of selected elements by ECM algorithm for geometric mapping φ AFFINE

Figure 20 :
Figure 20: Singular values decay profile of matrices of projected residuals, and number of selected elements by ECM algorithm for geometric mapping φ FFD+RBF

FOMFOMFOMFOMFigure 21 :Figure 22 :
Figure 21: QoI phase space plot for the ROMs against various HROMs for the training trajectory (Trajectory 1) for mapping φAFFINE

FOMFOMFOMFOMFigure 23 :
Figure 23: QoI phase plot for the ROMs against various HROMs for the training trajectory (Trajectory 1) for mapping

Fig. 24 Figure 24 :
Fig. 24 quantifies the trends observed in the phase space plots, by displaying the percentage errors incurred by the HROMs in the complete solution field and the QoI, when compared to the FOM and ROMs.

FOMFOMFOMFOMFigure 25 :Figure 26 :
Figure 25: QoI phase space plot for the ROMs against various HROMs for the testing trajectory (Trajectory 2) for mapping φAFFINE

FOMFOMFOMFOMFigure 27 :
Figure 27: QoI phase space plot for the ROMs against various HROMs for the testing trajectory (Trajectory 2) for mapping φFFD+RBF

Figure 28 :
Figure 28: Percentage error on QoI and solution field of ROM and FOM against various HROMs for the testing trajectory (Trajectory 2) for mapping φFFD+RBF

Figure 29 :
Figure 29: Speedup factors for the ROMs and HROMs presented in this example.The first column displays the performance of the ROMs, while the subsequent columns depict the HROMs' performance for each combination of truncation tolerances ϵSOL and ϵRES considered

Figure 30 :
Figure 30: Singular values of solution snapshots and the corresponding number of POD modes for both geometric mappings under Trajectory 1 and Trajectory 2 combined

Figure 31 :
Figure 31: QoI phase space plot for the FOM against various ROMs for the first training trajectory (Trajectory 2) for both geometric mappings

Figure 32 :
Figure 32: Percentage error on QoI and solution field of FOM against various ROMs for the first training trajectory (Trajectory 2) for both geometric mappings

Figure 33 :
Figure 33: QoI phase space plot for the FOM against various ROMs for the second training trajectory (Trajectory 1) for both geometric mappings

Figure 34 :
Figure 34: Percentage error on QoI and solution field of FOM against various ROMs for the second training trajectory (Trajectory 1) for both geometric mappings SOL =1e-03 SOL =1e-04 SOL =1e-05 SOL =1e-06 (a) Singular values for matrices Sr (ϵSOL) Number of selected elements for all combinations of ϵSOL and ϵRES

Figure 35 :Figure 36 :FOMFOMFOMFOMFigure 37 :
Figure 35: Singular values decay profile of matrices of projected residuals, and number of selected elements by ECM algorithm for geometric mapping φ AFFINE

Figure 38 :FOMFOMFOMFOMFigure 39 :
Figure 38: Percentage error on QoI and solution field of ROM and FOM against various HROMs for the first training trajectory (Trajectory 2) for mapping φAFFINE

Figure 40 :FOMFOMFOMFOM
Figure 40: Percentage error on QoI and solution field of ROM and FOM against various HROMs for the first training trajectory (Trajectory 2) for mapping φFFD+RBF

Figure 42 :
Figure 42: Percentage error on QoI and solution field of ROM and FOM against various HROMs for the second training trajectory (Trajectory 1) for mapping φAFFINE

FOMFOMFOMFOMFigure 43 :
Figure 43: QoI phase space plot for the ROMs against various HROMs for the second training trajectory (Trajectory 1) for mapping φFFD+RBF

Figure 44 :
Figure 44: Percentage error on QoI and solution field of ROM and FOM against various HROMs for the second training trajectory (Trajectory 1) for mapping φFFD+RBF

Figure 45 :
Figure 45: QoI phase space plot for the FOM against various ROMs for the testing trajectory (Trajectory 3) for both geometric mappings

Figure 46 :
Figure 46: Percentage on QoI and solution field of FOM against various ROMs for the testing trajectory (Trajectory 3) for both geometric mappings

FOMFOMFOMFOMFigure 47 :
Figure 47: QoI phase plot for the ROMs against various HROMs for the testing trajectory (Trajectory 3) for mapping φAFFINE

Figure 48 :
Figure 48: Percentage error on QoI and solution field of ROM and FOM against various HROMs for the testing trajectory (Trajectory 3) for φAFFINE

FOMFOMFOMFOM
Figure Percentage error on QoI and solution field of ROM and FOM against various HROMs for the testing trajectory (Trajectory 3) for mapping φFFD+RBF

Figure 50 :
Figure 50: Percentage error on QoI and solution field of ROM and FOM against various HROMs for the testing trajectory (Trajectory 3) for mapping φFFD+RBF

Figure 51 :
Figure 51: Speedup factors for the ROMs and HROMs presented in this example.The first column displays the performance of the ROMs, while the subsequent columns depict the HROMs' performance for each combination of truncation tolerances ϵSOL and ϵRES considered

Table 1 :
Type of mapping Number of Nodes Number of Elements Information of the mesh employed in this study for each geometric mapping.