LOR preconditioner for IncNS solver Pressure preconditioning

Issue/feature addressed

Please provide a description of the problem your changes address. Be sure to link to any related items in the issue tracker.

Improve on the existing pressure Poisson preconditioning available for the IncNS solver. For large 3D cases, the current default is Diagonal, which has high iteration counts with large standard deviations. Other preconditioners are not performant enough or do not scale well enough.

Proposed solution

A summary of the proposed changes and how they address the issue/feature at hand. Whenever possible, describe alternatives your considered and decisions you made.

Lower-order refined (LOR) preconditioner, also known as the SEMFEM preconditioner, uses a spectrally equivalent (Canuto et al. 2007) lower-order (P=1) discretisation to precondition the high-order problem - see figure below for an intuitive understanding. The motivation behind this approach is to achieve a sparser preconditioner by leveraging the sparsity of lower-order operators. It is also less memory-intensive and has minimal sensitivity to high aspect-ratio elements, making it suitable for complex geometries. A bounded iterative condition number behaviour with increasing problem size is also observed, meaning the larger the case, the more efficient this preconditioner is. This has been previously implemented for quads and hexes in libraries like Nek5000 and MFEM (Bello-Maldonado and Fischer 2019, Pazner 2020), but this is the first instance of implementation for prism and tetrahedrons. This is also the first time the LOR preconditioner is implemented for a modal expansion basis.

Implementation

A more detailed description of the changes. Please include any details necessary for reviewers to understand your work, and point out particular points would like them to pay particular attention to.

The following steps are done for the preconditioning:

• The input to the LOR preconditioner code is the residual of the iterative solution \hat{\mathbf{r}} at a given iteration (like all other preconditioners in Nektar++).
• The LOR space is generated at the PreconditionerLOR::v_BuildPreconditioner() step, using the MeshGraph::LinearMeshGraph method. The method has further methods setup to deal with the LOR space for each element type (LinMeshSetUp1DGeom, LinMeshSetUp2DGeom, LinMeshSetUpTetGeom, LinMeshSetUpPrismGeom). A coeffmap is also generated in this step, to deal with transformation from HO to LOR and vice versa (m_ho2lor).
• Interpolation of the solution from the HO discretisation to the LOR space is done using an appropriately constructed generalised Vandemonde matrix \mathcal{V}, depending on the choice of the higher-order expansion basis (MODIFIED/GLL_LAGRANGE). This is handled by MultiplyByBlockMatrix operations in the PreconditionerLOR::v_DoPreconditioner routine. Special keys (eGLLToCoeffs,eEquiSpacedToCoeffs) and their implementation have been set in the ExpList and StdExpansion to deal with this interpolation.
• After interpolating the LOR space, the linear system is solved using an appropriately chosen solution strategy, the inner-iteration strategy. For large 3D cases, the BoomerAMG solver is used via the PETSc interface. This is done using m_LORLinSys->SolveLinearSystem in PreconditionerLOR::v_DoPreconditioner.
• The solution from the LOR mesh is then copied back to the HO space, which is then transformed back onto the original expansion by the inverse of the \mathcal{V} matrix operation. Essentially, it is the inverse of the steps before the LOR system is solved.
• Appropriate scatter-gather transformations are handled using the assembly maps (m_ho2lor) in both spaces.

Tests

1. Helmholtz: Tests for LOR working for a Poisson solve in MultiRegions_Helmholtz3D with different element shapes: Triangles (Tri8_PreconditionerLOR_P8), Quads (Quad16_PreconditionerLOR_P6), Tetrahedrons (Tet_PreconditionerLOR_P8) and Prisms (Prism_PreconditionerLOR_P6).
2. IncNS solver: LOR working for both pressure and velocity systems in IncNavierStokesSolver for 3D elements only (Tetrahedrons, Prisms and Mixed meshes). Checked for XML and HDF5 formats. (Prism_channel_PreconditionerLOR, Tet_Kovasnay_PreconditionerLOR_pressure, Tet_Kovasnay_PreconditionerLOR_all, Tet_prism_channel_PreconditionerLOR_hdf5)
3. (Not completed) - Tests for MPI and using PETSc - need advice on how to proceed.

Suggested reviewers

Please suggest any people who would be appropriate to review your code. @ssherw, @dmoxey, @ccantwel

Notes

Please add any other information that could be useful for reviewers.

• To run this, IterativeFull and GMRES are needed.
• Dependencies: NEKTAR_USE_PETSC recommended when running for large 3D cases.
• LinMeshGraph routine is found to be memory intensive for large cases, purely due to the cost of constructing two different MeshGraphs
• PETSc options for using the BoomerAMG library need to be passed through the .petscrc file as of now
• For NUMMODES > 10, the electrostatic points distribution for the LOR Space in Triangles and Tetrahedrons starts causing errors, which causes the code to fail.
• To ensure proper prism ordering before LOR splitting in the LinMeshGraph routines, apply the ReorderPrisms module in NekMesh to the original mesh run with the LOR preconditioner. Warning: the ReorderPrisms module currently doesn't take care of the face curvature in the original mesh and removes that information. This becomes an issue with coarse meshes and high-curvature geometries.

Checklist

• Functions and classes, or changes to them, are documented.
• User guide/documentation is updated.
• Changelog is updated.
• Suitable tests added for new functionality.
• Contributed code is correctly formatted. (See the contributing guidelines).
• License added to any new files.
• No extraneous files have been added (e.g. compiler output or test data files).