ABL case files

group abl

Contains: user specified routines for the ABL case.

  • uservp() : variable properties

  • userf() : local acceleration term for fluid

  • userq() : local source term for scalars

  • userbc() : boundary conditions

  • useric() : initial conditions

  • userchk() : general purpose routine for checking errors etc.

  • userqtl() : thermal divergence for lowMach number flows

  • usrdat() : modify element vertices

  • usrdat2() : modify mesh coordinates

  • usrdat3() : general purpose routine for initialization

Functions

subroutine uservp(ix ix, iy iy, iz iz, ieg ieg)

Set variable properties, does not call any subroutines / functions.

Todo:

Implement turbulent Prandtl number parameter

param ix

x-index

param iy

y-index

param iz

z-index

param ieg

element index

subroutine userf(ix ix, iy iy, iz iz, eg eg)

Compute Coriolis acceleration.

\[ F_x, F_z = - f\hat{e}_y \times \vec{u} \]
See https://lists.mcs.anl.gov/pipermail/nek5000-users/2014-May/002798.html

Note

The Coroiolis acceleration also absorbs the part of the pressure under hydrostatic balance, which in turn drives the flow. When it is turned off, usrdat3() drives the flow.

param ix, iy, iz

[in] GLL point index

param eg

[in] global element number

subroutine userq(ix ix, iy iy, iz iz, eg eg)
subroutine useric(ix ix, iy iy, iz iz, eg eg)

Set up initial conditions.

subroutine userbc(ix ix, iy iy, iz iz, iside iside, eg eg)

Compute the boundary condition. See subdirectory bc/.

Note

it is OK to define both flux and temp here since only one would be used at the same time. See BCNEUSC(), BCDIRSC() and FACEIS()

subroutine userchk()

Compute the turbulent stress tensors and write statistics.

subroutine userqtl()
subroutine usrdat()
subroutine usrdat2()
subroutine usrdat3()

Compute inflow / outflow conditions a.k.a. driving force. Ubar=1 Not required when coriolis force is specified in userf.

See

https://github.com/Nek5000/NekExamples/blob/70a5792b04b7a4c2da16463f517863b10627398f/turbChannel/turbChannel.usr#L375-L386

Boundary conditions

group bc_moeng

Functions

subroutine abl_userbc(ix ix, iy iy, iz iz, iside iside, eg eg)

Stress boundary condition as formulated by Moeng (1984) Also has optional temporal filtering c.f: Yang et al. Physical Review Fluids 2, no. 10 (2017): 104601.

See

https://journals.ametsoc.org/doi/abs/10.1175/1520-0469(1984)041%3C2052:ALESMF%3E2.0.CO%3B2

See

https://doi.org/10.1103/PhysRevFluids.2.104601.

Note

Boundary condition evaluated at higher nodes, may need to be corrected for Ekmann turning

Note

This subroutine MAY NOT be called by every process

real function friction_vel(uh uh, kappa kappa, delta_z delta_z, z0 z0, Psi_M Psi_M)

Compute friction velocity.

real function thermal_flux(u_star u_star, kappa kappa, T_surf T_surf, T T, delta_z delta_z, z_os z_os, Psi_H Psi_H)

Complute thermal flux.

See

Eq. 14 in Gadde et al. (2020) doi:10.1007/s10546-020-00570-5

real function thermal_flux_fixed(x x)

Set thermal flux which varies along x, but is constant in time.

real function temp_strat(t_surf t_surf, y y, half_ymax half_ymax)

Set temperature dirichlet boundary conditions at bottom and top boundaries.

real function abl_pen_k(y_coord y_coord, z0 z0)

Compute the K term in penalty forcing.

param y_coord

mesh coordinate where the coefficient is evaluated

param z0

roughness length

group bc_channel

Variant of Moeng boundary condition which also allows to have a wall function at the top boundary.

Functions

subroutine abl_userbc(ix ix, iy iy, iz iz, iside iside, eg eg)

Stress boundary condition as formulated by Moeng (1984) Also has optional temporal filtering c.f: Yang et al. Physical Review Fluids 2, no. 10 (2017): 104601. https://doi.org/10.1103/PhysRevFluids.2.104601.

Note

Boundary condition evaluated at higher nodes, may need to be corrected for Ekmann turning

Note

This subroutine MAY NOT be called by every process

real function abl_pen_k(y_coord y_coord, z0 z0)

Compute the K term in penalty forcing.

param y_coord

mesh coordinate where the coefficient is evaluated

param z0

roughness length

group bc_noslip

Functions

subroutine abl_userbc(ix ix, iy iy, iz iz, iside iside, eg eg)

No slip and no penetration condition.

real function abl_pen_k(y_coord y_coord, z0 z0)

Compute the K term in penalty forcing.

Note

Be careful in specifying the penalty region, should be away from the wall

param y_coord

mesh coordinate where the coefficient is evaluated

param z0

roughness length

Forcing

group penalty_mini

A BC which appears as a/forcing”: the mini version.

Functions

subroutine pen_update()

Update penalty.

subroutine pen_forcing(ix ix, iy iy, iz iz, ieg ieg)

Compute penalty forcing.

Todo:

add ffy component and rotation

param ffx, ffy, ffz

[inout] forcing; x,y,z component

param ix, iy, iz

[in] GLL point index

param ieg

[in] global element number

subroutine pen_fmask_get()

Get 1D projection, array mapping and forcing smoothing.

This routine is just a simple version supporting only lines paralles to z axis. In future it can be generalised. The subroutine initializes pen_prj, pen_map and pen_npoint

See

Schlatter and Örlü, “Turbulent Boundary Layers at Moderate Reynolds Numbers.” pg. 12

Remark

This routine uses global scratch space CTMP0 and CTMP1

subroutine pen_fterms_get(ifreset ifreset)

Generate forcing terms: pen_k_len and pen_fstab.

This used to also initializes the pen_famp array. This was to allow for spatio-temporally varying forcing, which has been removed in this version.

param ifreset

[in] reset flag to reinitialize geometry term pen_k_len

subroutine pen_register()

Register penalty module.

Note

This routine should be called in frame_usr_register

subroutine pen_init()

Initilise penalty module.

Note

This routine should be called in frame_usr_init

logical function pen_is_initialised()

Check if module was initialised.

return

pen_is_initialised

subroutine local_cart_dy(us us, u u, N N, e e, Dt Dt)

Assuming the mesh is cartesian, compute exact y-derivative for an element, by operating.

\[ \mathbf{u} D^T \]

param us

[out]

param u, N, e, Dt

[in]

subroutine local_cart_dy_pseudo(us us, u u, N N, e e, Dt Dt)

Assuming the mesh is cartesian, compute pseudo(?) y-derivative along y-axis for an element, by operating.

\[ Dt \mathbf{u} \]

param us

[out]

param u, N, e, Dt

[in]

subroutine cart_diff_y(uy uy, u u, ifexact ifexact)