Implicit Incompressible Solver using Velocity-Correction scheme

Issue/feature addressed

Simulations of convection-dominated flows (high Reynolds number) with the Semi-Implicit VelocityCorrectionScheme (VCS) often have strong restrictions on the time step size. This restrictions requires users to set significantly smaller time-step sizes due to stability rather than accuracy of the physical solution. Consequntly, the user must deal with higher computational cost to reach large physical times.

Proposed solution

This MR adds an implicit scheme for the IncNavierStokesSolver to Nektar++. This scheme allows the user to choose larger time-step sizes in exchange for a cost-overhead. It is based on the work by Dong and Shen (2010).

Implementation

The implicit scheme VCSImplicit is derived from the VelocityCorrectionScheme class. It uses most of the existing functionality for calls from the TimeIntegrationAlgorithmGLM. The time integration calls DoOdeRhs() which corresponds to EvaluateAdvection_SetPressureBcs(), the explicit part of the VCS, and DoImplicitSolve() which corresponds to SolveUnsteadyStokes(). SolveUnsteadyStokes() computes forcing terms as well as the solves for pressure and velocity.

Compared to the original VCS, VCSImplicit has the main difference that it does not require extrapolation of the explicit advection terms and uses the weak pressure form that is implemented in VCSWeakPressure. So, it reimplements EvaluateAdvection_SetPressureBcs() and evaluates the advection terms without handing them back to the TimeIntegrationAlgorithm for extrapolation (via outarray). This causes some unused memory which is not a problem with the current smaller scale tests, but might need a better solution later. The subsequent call to EvaluatePressureBcs is similarly modified to not use extrapolated advection and copy functionality from WeakPressureExtrapolation. To disable extrapolation, it was necessary to derive a new class ImplicitExtrapolation from WeakPressureExtrapolation and remove the extrapolation in EvaluatePressureBcs.

Moving to the SolveUnsteadyStokes routine, the main difference to VelocityCorrectionScheme is again removing the extrapolated advection and body forcing term usually included in inarray as: \sum_q \alpha_q u^{n-q} + dt*(f^{n+1} - \sum_q \beta_q N(u)^{n-q}). Consequently, we have to manually add f^{n+1} and the advection N(u)^n to the pressure and velocity forcing terms in SetupPressureForcing() and SetupViscousForcing(), respectively. Additionally, the advecction velocity \tilde{u} is computed within SetupViscousForing(). In general, the pressure routines closely follow the VCSWeakPressure implementation. For the velocity solve, the SolveViscous() call is modified to, first, add the advection velocity as a variable coefficient (varcoeff) and, secondly, use a LinearAdvectionDiffusionReactionSolve() instead of the HelmSolve(). Due to conventions in the subsequent routines, the sign of the lambda coefficient is flipped.

In comparison to the original VCS, the implicit solve uses a matrix with time-varying coefficients. Thus, the global matrix system needs to be deleted after each time step. This is done using the UnsetGlobalLinSys() routine which is called after each new velocity has been solved for. This routine is strictly necessary to avoid a memory leak.

Minor changes allow non-symmetric Laplacian matrices via varcoeffs. This branch implements a switch that detects matrix keys of type: AdvectionDiffusionReaction, throws a warning and changes the iterative solver to GMRES, if it was not already selected.

Tests

• Test for stability of the scheme on 2D Kovasznay flow with 4-Quad mesh using polynomial order P=9 and P=15
• Check for Taylor-Hood expansions on 2D Kovasznay flow
• Tests for the HOutflow and convective-like boundary conditions (Based on tests for VelocityCorrectionScheme)

Suggested reviewers

• David Moxey
• Spencer Sherwin

Notes

• The scheme is currently implemented and mostly tested for 2D.
• There was an issue with making SolveUnsteadyStokes virtual. The SmoothedProfileMethod is derived from VelocityCorrectionScheme and constructs a virtual function of SolveUnsteadyStokes. This merge resolves this by defining a virtual SolveUnsteadyStokes within VelocityCorrectionScheme instead, see commit 61f1b6f9
• There was an issue with the GMRES-Switch and the LinearElasticSolver test of L-Domain.xml. The switch now detects only AdvectionDiffusionReaction, instead of any Advection type matrix, because the LinearElasticSolver uses an AdvectionReaction key, but the L-Domain test fails when it does not run on a CG solver. There is likely a better method to detect a non-symmetric matrix in Nektar (e.g. the MatrixStorageType?)
• There is room for improvement, however, this first merge should allow the solver to evolve more inline with current developments on the master branch, instead of working on a year-old branch.

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).