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 forVelocityCorrectionScheme
)
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. TheSmoothedProfileMethod
is derived fromVelocityCorrectionScheme
and constructs a virtual function ofSolveUnsteadyStokes
. This merge resolves this by defining a virtualSolveUnsteadyStokes
withinVelocityCorrectionScheme
instead, see commit 61f1b6f9 - There was an issue with the GMRES-Switch and the
LinearElasticSolver
test ofL-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. theMatrixStorageType
?) - 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).