In this tutorial we re-visit the solution of fluid-structure interaction problems with (pseudo-)solid fluid mesh updates, focusing on the efficient solution of the governing equations by GMRES, using the problem-specific block-preconditioner developed in
where full details of the analysis and further results from a variety of other test problems can be found.
We will demonstrate the development and application of the preconditioner using the problem of flow in a 2D channel with an elastic leaflet discussed in another tutorial and illustrated in the sketch below:
Flow is driven (via an imposed parabolic inflow profile) through a 2D channel which is partially obstructed by a thin-walled elastic leaflet. The fluid leaves the domain at the far downstream end of the channel where we assume parallel, axially traction-free outflow.
We use the same discretisation (2D quadrilateral Taylor-Hood elements for the fluid; 1D Hermite beam elements for the leaflet) as in the original tutorial. However, here we perform the fluid mesh update in response to the deformation of the elastic leaflet using a pseudo-elastic approach. This is done by "wrapping" the underlying
QTaylorHoodElement into a
PseudoSolidNodeUpdateElement and employing
ImposeDisplacementByLagrangeMultiplierElements to impose the deformation of the elastic leaflet onto the (pseudo-solid) fluid mesh. The following tutorials discuss the relevant methodologies:
The discretised problem contains the following types of discrete unknowns:
Using this classification of the unknowns, the linear system to be solved in the course of the Newton iteration can be re-ordered into the following block structure:
where the vector contains the Navier-Stokes (fluid) unknowns (velocities and pressure), represents unknowns describing the deformation of the fluid-loaded solid, represents the nodal positions in the fluid mesh, and contains the discrete Lagrange multipliers which impose the deformation of the FSI boundary in the fluid mesh. The three diagonal blocks in the Jacobian matrix, , are the two ``single-physics'' Jacobian matrices ( , the Navier-Stokes Jacobian; , the tangent stiffness matrix of the fluid-loaded solid) and the Jacobian,
associated with the Lagrange-multiplier-constrained pseudo-elasticity problem governing the fluid mesh update, discussed in more detail in another tutorial. The non-zero off-diagonal blocks arise through the interactions between the various single-physics problems: represents the effect of the pseudo-solid unknowns (the nodal positions in the fluid mesh) on the discretised Navier–Stokes equations – this incorporates the so-called shape-derivatives and the time-discretised version of the no-slip condition on the FSI boundary which expresses the fluid velocities in terms of the nodal velocities of the fluid mesh. captures the effect of the fluid traction (pressure and shear stresses) on the fluid-loaded solid. Since the shear stresses depend on gradients of the fluid velocity, the traction is also affected by the nodal positions in the fluid mesh – this dependency gives rise to the matrix . Finally, arises because the prescribed boundary displacement of the fluid mesh depends on the displacement field in the fluid-loaded solid. The zero block reflects the fact that the discretised Navier-Stokes equations do not depend directly on the displacement field of the actual solid. The zero block indicates that the pseudo-solid equations are not affected by the fluid unknowns. The block is zero because the real solid affects the pseudo-solid only indirectly via the Lagrange multiplier constraint (which gives rise to ). Finally, the non-zero blocks in the last block column reflect the fact that the Lagrange multiplier is only used to enforce the displacement constraint in a one-way coupling in which the real solid affects the pseudo-solid but not vice versa. See Muddle, Mihajlovic & Heil (2012) for more detail.
To construct a preconditioner for the linear system (1) we replace the bottom-right block in the Jacobian with the pseudo-solid preconditioner,
which we discussed in another tutorial. This yields the pseudo-solid FSI preconditioner
The preconditioner can be shown to have a block triangular structure. Its application therefore requires the solution of four subsidiary linear systems involving the matrices , , and and four sparse matrix-vector products with the interaction matrices , , and .
Numerical experiments in Muddle, Mihajlovic & Heil (2012) show that an efficient implementation of the preconditioner is obtained by
and by replacing the remaining block-solves within these preconditioners by a small number of AMG cycles or CG iterations. We retain SuperLU as the (exact) direct solver for the linear system involving the real solid's tangent stiffness matrix .
With these approximations, the computational cost of one application of the preconditioner is linear in the number of unknowns (see Comments and Exercises for a more detailed discussion of this issue). The optimality of the preconditioner can therefore be assessed by demonstrating that the number of GMRES iterations remains constant under mesh refinement.
The preconditioner described above is implemented within
oomph-lib's (parallel) block preconditioning framework which is described in another tutorial. For the purpose of the implementation, we decompose the preconditioning matrix into the 3x3 main blocks indicated by the vertical and horizontal lines in (1) and (3).
The degrees of freedom within these sub-blocks are classified into "dof-types" on an element-by-element basis, using the functions
GeneralisedElement::ndof_types(...) which returns the number of dof types classified by that element.
GeneralisedElement::get_block_numbers_for_unknowns(...) which associates each degree of freedom (identified by its global equation number) with the dof type within its block.
The default implementation of these functions within the Navier-Stokes elements (which differentiate between the fluid velocity components and the pressure) and the beam elements (which do not distinguish between the different types of degree of freedom) is appropriate for the use with the (subsidiary) preconditioners employed here. The degrees of freedom representing the nodal positions in the (pseudo-solid) fluid mesh need to be sub-divided into those that are and are not constrained by the Lagrange multipliers that impose the deformation of the FSI boundary. This is most conveniently done by overloading the relevant functions in a templated wrapper element,
PseudoElasticBulkElement<...>, as discussed in another tutorial.
We examine the performance of the preconditioner in steady and unsteady test problems. First we perform a sequence of steady solves, incrementing the Reynolds number in steps of 25. The final steady solution is then used as the initial condition for an unsteady simulation in which the inflow is subjected to a time-periodic oscillation. The tables below show the GMRES iteration counts (averaged over all linear solves performed in the course of all Newton iterations) as a function of the mesh refinement (represented by the total number of unknowns) for different implementations of the preconditioner.
|GMRES (blocks solved by SuperLU)||5.88889||5.88889||5.88889||6.11111|
|GMRES (LSC; blocks solved by SuperLU)||27.8889||36.1111||42.3333||47.2222|
|GMRES (LSC; pseudo-solid; blocks solved by Hypre/CG)||37.3333||43.2222||49.3333||55.3333|
|GMRES (blocks solved by SuperLU)||6.5124||6.95868||7.15323||7.2619|
|GMRES (LSC; blocks solved by SuperLU)||15||15.9256||16.2097||16.6429|
|GMRES (LSC; pseudo-solid; blocks solved by Hypre/CG)||22.2066||25.3967||27.0484||27.9921|
As expected from the theory (see Muddle, Mihajlovic & Heil (2012)), the GMRES iteration counts are small and virtually mesh independent for the exact implementation of the preconditioner. Replacing the (costly) exact solves by (faster) approximate solves leads to a modest increase in the (absolute) number of GMRES iterations. The most significant deterioration of the iteration counts results from the use of the LSC Navier Stokes preconditioner as the inexact solver (subsidiary preconditioner) for the Navier-Stokes block. The behaviour observed here (slight mesh dependence for steady solves; virtual mesh independence (once the mesh is sufficiently fine) for unsteady solves) mirrors that observed in single-physics Navier-Stokes problems.
The benefit of switching to the approximate solves becomes apparent in the next two tables which shows the average cpu times required for the solution of the linear systems by GMRES. (For comparison the times labelled as "SuperLU" show the times when the linear systems (1) are solved directly by SuperLU rather than by GMRES; they illustrate how quickly direct solvers become uncompetitive.) For sufficiently fine discretisations the larger number of GMRES iterations for the inexact solves is more than compensated for by the much lower cost of the preconditioning operations. The final implementation yields cpu times that (at least for unsteady problems) are proportional to the number of unknowns – the hallmark of an optimal solver.
|GMRES (blocks solved by SuperLU)||2.31567||15.1042||58.5491||147.531|
|GMRES (LSC; blocks solved by SuperLU)||2.8051||18.3197||58.7335||138.558|
|GMRES (LSC; pseudo-solid; blocks solved by Hypre/CG)||3.47733||15.1812||37.6856||73.0526|
|GMRES (blocks solved by SuperLU)||1.73668||11.1679||42.2051||105.97|
|GMRES (LSC; blocks solved by SuperLU)||1.87403||11.0113||34.2697||81.4382|
|GMRES (LSC; pseudo-solid; blocks solved by Hypre/CG)||2.26537||9.40194||22.5644||41.8251|
Most of the driver code is unchanged from the implementation discussed in the original tutorial. We therefore only discuss the modifications required
To perform the node update in the fluid mesh in response to the deformation of the leaflet using the equations of large-displacement elasticity we "wrap" the Navier-Stokes element in the doubly-templated
PseudoSolidNodeUpdateElement<WRAPPED_ELEMENT, SOLID_ELEMENT> class. The first template argument of this class,
WRAPPED_ELEMENT, specifies the type of the "wrapped" element – here the nine-node, 2D quadrilateral
QTaylorHoodElement<2>. The second template argument,
SOLID_ELEMENT, specifies the
SolidElement used to perform the node update – here the nine-node, 2D quadrilateral
QPVDElement<2,3>. Since the preconditioner requires the (re-)classification of the pseudo-solid degrees of freedom into constrained and unconstrained nodal positions, we wrap this element too, using the
PseudoElasticBulkElement already discussed in another tutorial. The full specification of the fluid element is therefore given by
The main code only requires minor changes, all associated with the specification of different solver options. Since we provide the option to use Hypre and Trilinos solvers, we need to activate MPI if
oomph-lib has been compiled with MPI support (even if the code is run in serial):
Next we define and process the possible command line flags which are used to select solver options and to specify the spatial resolution
We then assign the command line flags that were specified and document the ones that were recognised, before building the problem, specifying the heavily wrapped version of the
Navier-Stokes element as the template parameter.
The rest of the main code is very similar to the original version discussed in another tutorial and is therefore omitted here.
We provide a new function,
set_iterative_solver() to switch the solver to GMRES and to instantiate the FSI preconditioner. We start by creating the GMRES solver, using Trilinos' implementation if it is available:
Next, we create an instance of the FSI preconditioner for a 2D problem and pass it to the solver:
The classification of the degrees of freedom within the block preconditioning framework requires the specification of the (pseudo-solid) fluid mesh, the "real" solid mesh and the mesh containing the Lagrange multiplier elements that impose the deformation of the FSI boundary in the fluid mesh:
We provide the option to bypass the use of the LSC Schur complement preconditioner as the subsidiary preconditioner (inexact solver) for the Navier-Stokes block:
By default the preconditioner performs the preconditioning operations as described above, using SuperLU for the block solves in (3). It is possible to specify alternative (approximate) solvers for these. For instance, the behaviour of the LSC preconditioner can be modified by following the pointer to the Schur complement preconditioner
whose behaviour can be modified as discussed in LSC preconditioner tutorial. For instance, to use a block-triangular approximation for the Navier-Stokes momentum block we specify
If Hypre is available, we can then (approximately) solve the diagonal blocks using a small number of AMG V-cycles (with settings specified in the helper function
Similarly, we can use Hypre (with settings appropriate for a 2D Poisson problem) to approximately solve the pressure Poisson systems in the computations of the Schur complement approximation:
Approximate solvers for the diagonal blocks in the preconditioner (2) for the the Lagrange-multiplier constrained pseudo-elastic equations can be chosen as discussed in another tutorial. For instance, we can employ a block-upper triangular approximation for the augmented pseudo-elastic block :
Its diagonal blocks can then again be solved by Hypre,
Similarly, the linear systems involving the mass matrices in can be solved using diagonally preconditioned CG:
run.shused to generate the data presented above, may be helpful.
A pdf version of this document is available.