Notes on Nek5000 internals#
Nek5000 core library#
Note
The following is auto-generated using Doxygen and represents Nek5000 v19. The official Nek5000 documentation is hand-written and is more general.
External links#
The FLOW-SESE Course 2021 book - Nek5000: Theory, Implementation and Optimization
Boundary condition#
Algorithm#
The following describes the flow of logic in setting the boundary conditions.
In drive1.f
:
nek_advance()
: single timestep if Pn/Pn-2igeom=1
: RHS termsigeom=2
: Actualcall heat
call fluid
-> drive2.f
fluid()
:iftran=1,ifrich=0
->plan3()
inplanx.f
plan3()
: Helmholtz operator split into two stepsigeom=1,2
makef()
: extrapolate pressure, solve for intermadiate velocitysolve for update of pressure, velocity divergence free. Only approximate operator is used here, because the exact one requires Helmholtz operator.
Block LU decomposition -> Consistent Helmholtz operator
Boundary condition is not required for the second step. Preconditioner.
See Blair Perot (1993)
Also if igeom=2
:
sethlm
:h1=diffusion
h2=helmholtz operator
cresvif()
compute residual:lagvel, save velocity field
bcdirvc
: boundary dirichlet velocity,bcneutr
: boundary neutral ->bdry.f
opgradt
: extrapolate pressire in timeadd residual
ophinv
,opadd2
: inverse helmholtzincomprn()
:Project velocity
opgradt
: pressureopbinv
->navier1
:bdry?
opmask
: set 0 to dirichlet BC, masked from the update
Boundary condition subroutines#
Periodicity: node numbering
bcdirvc
: Dirichlet#
CB = V (Constant Dirichlet in rea file)
CB = v (user provided Dirichlet in userbc)
:FACEIV
: works on temporary arraystmp1, tmpr3
facind: face index
facev
->userbc
opdsop
: max min, to get consistent bc shared nodesif not stress forulation
ifstrs
:add2
-> Velocity + Tmp
bcneutr
: Neumann transient?#
No mask
chkcbc
: alignment checkSYM
<- calledbcmask
<->=nek_init
modifies the masking for one componentAlso by
sethmsk
b1mask
,b2mask
,b3mask
arrays:FEM trick, if you don’t anything: Natural BC. Normal stresses are 0, BC is the fluxes on the last element.
subs2.f
sh
:faceiv
-> sets traction array TRX, TRY, TRZ -> as accelerationfaccvs
-> multiplies by the surface area -> forceglobrot
, does rotation with respect to global coordinatesadded to RHS:
bfx
,bfy
-> acceleration
Velocity should be v=0
Caution: Pressure preconditioner and BC are related
Todo
Set velocity v=0, normal component of velocity is zero, opposite of on
boundary condition.
Time-stepping algebra#
Mass matrix requires time levels (n, n-1, n-2): extrapolated Conv / Stiffness requires time levels (n-1, n-2): explicit
Comparison with Adams-Bashforth#
AB2: needs previous RHS BDF: needs previous LHS
Why not: Runge-Kutta?#
According to mathematicians, it is not clear because you have a pressure correction in each sub-step.
CFL number#
BDF3/EXT3: CFL=0.6. Characteristic method
Filtering#
Filtering in Nek5000 is controlled by
We prefer to use HPFRT (high-pass filtering with relaxation term) in Nek5000. There is also an
HPFRT vs explicit filtering#
Similarities#
Both methods aims to extract energy from higher modes of the solution.
Filtering weight determines the strength of the filter
Filter cut-off ratio determines the what part of the spectrum is diffused numerically by the filter,
Differences#
HPFRT introduces a term in the RHS of the governing equations, thereby not violating zero velocity divergence, or in other words behaving like a source or sink.
Explicit filtering is performed at the end of the time step, and very likely could result in velocity divergence.
Filtering as a sponge-layer#
Filter weight can be varied within the domain, and this would allow the filter to behave like a sponge-layer if needed.