exabl / snek5000-abl

Nek5000 ABL user code

snek5000-abl

Simulations of turbulent atmospheric boundary layer using snek5000.

Documentation: https://exabl.github.io/snek5000-abl/

Warning

:::{warning} The code is a prototype and far from ready for production runs. Some parts of the code are well validated and some parts require rigorous testing and could contain bugs. Use it with caution. To know the detailed status of the code, checkout the roadmap. :::

Quick start

Install using Python 3.8+ as follows:

git clone --recursive https://github.com/exabl/snek5000-abl.git
cd snek5000-abl
pip install -e .

Activate necessary environment variables

source activate.sh

Use the command line tool to launch / set / inspect the simulation parameters

abl --help

Tests

pip install -e '.[tests]'
# Run simple tests: including compilation
pytest
# Run slow tests: launches simulation
pytest --runslow

snek5000-abl's Issues

Measuring Obukhov length

Suggestion from @vezeli

Here is an idea for you if you have the time and interest to follow it... When you have your LES result you could look at G/P (ratio of buoyancy production and shear production) this should give value of Obukhov length scale, than you can compare it with the value obtained from the standard definition of Obukhov length. The difference is quite substantial because in one case you estimate the Obukhov length from surface scaling properties while in the other formulation you use turbulence fluxes throughout the atmospheric boundary layer 🤓

This can be checked when we do stable ABL simulations.

Warnings while restarting and quickly goes NaN

While using revision 1054a63

WARNING: restart file has a NSPCAL > LDIMT
read only part of the fld-data!
WARNING: NPSCAL read from restart file differs from
currently used NPSCAL!

Alignment in common block

Ensure memory alignment by using same datatypes in user defined common blocks

Possible improvements to SGS implementation

• Average $nu_t$ within the element
• Compute $nu_t$ and dU from the log-law. And modify viscosity instead of setting shear stress at the wall (Timofey's approach)
• Look at spectra (requires #41) and test decaying field turbulence. Is it k-5/3 till grid size and cutoff? There should not be any energy pile up

Postprocessing utilities in the package

• spatial_means: pandas load
• stats: pymech load, plot different statistics
• phys_fields: pymech load, shenfun operator, spectra computation.

Explore different SGS models

• Dynamic Smagorinsky Model: present in turbChannel example in NekExamples. See also: Nek5000/NekExamples#105. More computationally expensive, may not be needed. Wall damping function is not hardcoded, and but automatic / dynamic

• Algebraic model for computing unresolved K.E + Deodroff.

• Vreman model?: no wall damping required.

• Piomelli Guertz: See https://doi.org/10.1017/jfm.2015.29 . Grid scale independent SGS model, perhaps useful for spectral element methods, because non-uniform grid.

• Shear improved Smagorinsky model: See https://doi.org/10.1017/S0022112006003429. Simple to implement, take away the mean shear component and only fluctuation shear is retained.

• Algebraic Minimum Dissipation (AMD) model: Rozema et al., “Minimum-Dissipation Models for Large-Eddy Simulation.”

Random seed

Find a way to uniquely set random seed for a given mpi rank and size, while not resulting in similar sequences

Buoyancy / stratification

• issue with turning on variableProperties
• top radiative B.C / sponge layer
• bottom Moeng B.C. with correction terms

Upgrade to v19

Also use new statistics toolbox for v19.

CI: Fix docs workflow

Exception occurred:
  File "/opt/hostedtoolcache/Python/3.7.7/x64/lib/python3.7/site-packages/breathe/parser/__init__.py", line 93, in parse
    raise ParserError(e, filename)
breathe.parser.ParserError: file /home/runner/work/snek5000-abl/snek5000-abl/docs/_build/xml/makefile__usr_8inc.xml: mismatched tag: line 61, column 33

Alternate stress boundary condition evaluation

Timofey suggested two differences on how it is done in Fluid mechanics:

• u* is computed using Newton Raphson method by evaluating u+, y+ curve
• Specify no slip boundary condition and manipulate $\nu$ at wall to mimic roughness because $\tau = \nu \Delta u$.

Setup debug flags and document it

Debugging should be capable of:

• stepping through the code and inspecting values of it
• flag to stop execution when any values hit NaN
• set up a breakpoint just before / after things blow up and goes NaN

Noise in initial condition

Following #30 (specifically at 4ebe893 ), for test purposes initial condition noise was commented out. Adjust it to moderate levels and put it back in.

Fix Doxygen output

Currently it fails when the src/abl/bc directory is included

Modelling improvements?

• clean freestream initial condition

• constant SGS filter length scale ~ maximum grid spacing, $\Delta$, decrease $C_0$ in Smagorinsky model to 0.15 (#21)

• Boundary condition: evaluate shear stress $\tau_w$ based on velocities at 3rd and 4th grid levels (#20)

• Grid refinement: Non uniform elements in the vertical direction (master)

• constant SGS viscosity sans wall damping

• SGS wall damping function: Change power from 0.5 to 3.0 (#22)

Should add parameters to toggle this?

Temperature-related fields in statistics

Hi @ashwinvis, hope things are going well! I was wondering if you've implemented adding temperature-related fields to the statistics.
I checked here
https://github.com/exabl/snek5000-abl/blob/master/src/abl/toolbox/stat_IO.f
and it seems like it's not there at least. Thanks!

Floating point exception

In branch enh/simple-postproc. cc: @geert-brethouwer

Traceback

Backtrace for this error:

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0  0x7f6e07821faf in ???
exabl/eturb#1  0x55aac369a75d in gather_double_add
        at src/gs_local.c:104
exabl/eturb#2  0x55aac369e277 in gslib_gs_gather
        at src/gs_local.c:211
exabl/eturb#3  0x55aac36911f1 in cr_exec
        at src/gs.c:612
exabl/eturb#4  0x55aac3692b7b in dry_run_time
        at src/gs.c:977
exabl/eturb#5  0x55aac3692e52 in auto_setup
        at src/gs.c:1018
exabl/eturb#6  0x55aac36934aa in gs_setup_aux
        at src/gs.c:1113
exabl/eturb#7  0x55aac3693992 in fgslib_gs_setup_pick_
        at src/gs.c:1202
exabl/eturb#8  0x55aac3693a07 in fgslib_gs_setup_
        at src/gs.c:1210
#0  0x7fd858147faf in ???
exabl/eturb#9  0x55aac3469a09 in setupds_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/dssum.f:20
exabl/eturb#1  0x5560850b575d in gather_double_add
        at src/gs_local.c:104
exabl/eturb#2  0x5560850b9277 in gslib_gs_gather
        at src/gs_local.c:211
#0  0x7fd4e8cbdfaf in ???
exabl/eturb#3  0x5560850ac1f1 in cr_exec
        at src/gs.c:612
exabl/eturb#4  0x5560850adb7b in dry_run_time
        at src/gs.c:977
exabl/eturb#5  0x5560850ade52 in auto_setup
        at src/gs.c:1018
exabl/eturb#6  0x5560850ae4aa in gs_setup_aux
        at src/gs.c:1113
exabl/eturb#7  0x5560850ae992 in fgslib_gs_setup_pick_
exabl/eturb#10  0x55aac345b893 in setup_topo_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/connect1.f:84
        at src/gs.c:1202
exabl/eturb#8  0x5560850aea07 in fgslib_gs_setup_
        at src/gs.c:1210
exabl/eturb#1  0x56206cb3d75d in gather_double_add
        at src/gs_local.c:104
exabl/eturb#2  0x56206cb41277 in gslib_gs_gather
        at src/gs_local.c:211
#0  0x7ff23fb6dfaf in ???
exabl/eturb#11  0x55aac3427e1b in nek_init_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive1.f:93
exabl/eturb#12  0x55aac34275ed in nekton
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive.f:12
exabl/eturb#13  0x55aac3427669 in main
exabl/eturb#3  0x56206cb341f1 in cr_exec
        at src/gs.c:612
exabl/eturb#4  0x56206cb35b7b in dry_run_time
        at src/gs.c:977
exabl/eturb#5  0x56206cb35e52 in auto_setup
        at src/gs.c:1018
exabl/eturb#6  0x56206cb364aa in gs_setup_aux
        at src/gs.c:1113
exabl/eturb#7  0x56206cb36992 in fgslib_gs_setup_pick_
        at src/gs.c:1202
exabl/eturb#8  0x56206cb36a07 in fgslib_gs_setup_
        at src/gs.c:1210
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive.f:18
exabl/eturb#1  0x55dba92fb75d in gather_double_add
        at src/gs_local.c:104
exabl/eturb#9  0x556084e84a09 in setupds_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/dssum.f:20
exabl/eturb#2  0x55dba92ff277 in gslib_gs_gather
        at src/gs_local.c:211
#0  0x7f0b02067faf in ???
exabl/eturb#9  0x56206c90ca09 in setupds_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/dssum.f:20
exabl/eturb#3  0x55dba92f21f1 in cr_exec
        at src/gs.c:612
exabl/eturb#4  0x55dba92f3b7b in dry_run_time
        at src/gs.c:977
exabl/eturb#5  0x55dba92f3e52 in auto_setup
        at src/gs.c:1018
exabl/eturb#1  0x55a36b6e575d in gather_double_add
exabl/eturb#6  0x55dba92f44aa in gs_setup_aux
        at src/gs.c:1113
exabl/eturb#7  0x55dba92f4992 in fgslib_gs_setup_pick_
        at src/gs.c:1202
exabl/eturb#8  0x55dba92f4a07 in fgslib_gs_setup_
        at src/gs.c:1210
        at src/gs_local.c:104
exabl/eturb#2  0x55a36b6e9277 in gslib_gs_gather
        at src/gs_local.c:211
exabl/eturb#10  0x556084e76893 in setup_topo_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/connect1.f:84
exabl/eturb#3  0x55a36b6dc1f1 in cr_exec
        at src/gs.c:612
exabl/eturb#4  0x55a36b6ddb7b in dry_run_time
        at src/gs.c:977
exabl/eturb#5  0x55a36b6dde52 in auto_setup
        at src/gs.c:1018
exabl/eturb#6  0x55a36b6de4aa in gs_setup_aux
        at src/gs.c:1113
exabl/eturb#7  0x55a36b6de992 in fgslib_gs_setup_pick_
        at src/gs.c:1202
exabl/eturb#8  0x55a36b6dea07 in fgslib_gs_setup_
        at src/gs.c:1210
exabl/eturb#11  0x556084e42e1b in nek_init_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive1.f:93
exabl/eturb#12  0x556084e425ed in nekton
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive.f:12
exabl/eturb#13  0x556084e42669 in main
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive.f:18
exabl/eturb#10  0x56206c8fe893 in setup_topo_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/connect1.f:84
exabl/eturb#9  0x55dba90caa09 in setupds_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/dssum.f:20
exabl/eturb#9  0x55a36b4b4a09 in setupds_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/dssum.f:20
exabl/eturb#11  0x56206c8cae1b in nek_init_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive1.f:93
exabl/eturb#12  0x56206c8ca5ed in nekton
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive.f:12
exabl/eturb#13  0x56206c8ca669 in main
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive.f:18
exabl/eturb#10  0x55dba90bc893 in setup_topo_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/connect1.f:84
exabl/eturb#11  0x55dba9088e1b in nek_init_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive1.f:93
exabl/eturb#12  0x55dba90885ed in nekton
exabl/eturb#10  0x55a36b4a6893 in setup_topo_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/connect1.f:84
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive.f:12
exabl/eturb#13  0x55dba9088669 in main
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive.f:18
exabl/eturb#11  0x55a36b472e1b in nek_init_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive1.f:93
exabl/eturb#12  0x55a36b4725ed in nekton
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive.f:12
exabl/eturb#13  0x55a36b472669 in main
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive.f:18
#0  0x7f353e0aefaf in ???
exabl/eturb#1  0x559c8b34875d in gather_double_add
        at src/gs_local.c:104
exabl/eturb#2  0x559c8b34c277 in gslib_gs_gather
        at src/gs_local.c:211
exabl/eturb#3  0x559c8b33f1f1 in cr_exec
        at src/gs.c:612
exabl/eturb#4  0x559c8b340b7b in dry_run_time
        at src/gs.c:977
exabl/eturb#5  0x559c8b340e52 in auto_setup
        at src/gs.c:1018
exabl/eturb#6  0x559c8b3414aa in gs_setup_aux
        at src/gs.c:1113
exabl/eturb#7  0x559c8b341992 in fgslib_gs_setup_pick_
        at src/gs.c:1202
exabl/eturb#8  0x559c8b341a07 in fgslib_gs_setup_
        at src/gs.c:1210
exabl/eturb#9  0x559c8b117a09 in setupds_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/dssum.f:20
exabl/eturb#10  0x559c8b109893 in setup_topo_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/connect1.f:84
exabl/eturb#11  0x559c8b0d5e1b in nek_init_
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive1.f:93
exabl/eturb#12  0x559c8b0d55ed in nekton
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive.f:12
exabl/eturb#13  0x559c8b0d5669 in main
        at /home/avmo/src/exabl/eturb/lib/Nek5000/core/drive.f:18

Fix planar average

According to Adam:

... After brief check. Your implementation is wrong, as you mess up element internal degrees of freedom with mesh structure. If I understand correctly routine you are trying to use is based on subroutine planar_avg(ua,u,hndl) in navier5.f. This one refers to subroutine gtpp_gs_setup(hndl,nelgx,nelgy,nelgz,idir) in dsssum.f However the meaning of the argument refers to the global mesh structure, not to the number of degrees of freedom inside single element. In this routine:
hndl - flag referring to proper gssum initialisation
nelgx,nelgy,nelgz - number of elements in x,y,z direction (similar to line parameter (lelx=1,lely=1,lelz=1) ! global tensor mesh dimensions; in SIZE file)
idir - averaging direction.

This routine has significant limitations as it can be used to meshes build with genbox. There can be exceptions, but gtpp_gs_setup can be used for box geometry with proper element numbering (like in genbox). It is not enough to have uniform direction. It has to be a box which is simple tensor product. In your code you write
call gtpp_gs_setup(igs_x, nx1, ny1, nz1, 1) ! x-avg
This is wrong, as nx1, ny1, nz1 give number of degrees of freedom inside element, not number of elements in each direction. This is the main reason routine does not work. You could define your mesh structure (if it is a box) in SIZE and later reference lelx, lely, lelz calling gtpp_gs_setup.

Required for #25

Change coordinate system to make z the wall normal direction

On the Nek level, it would need changes in:

• box file
• various places in abl.usr and smagorinsky.f

On the snek level, only the parameters need to be changed?

Boundary condition for SGS models

As shown in Eqns 19 & 20 in the paper

Gadde, Stieren, and Stevens, “Large-Eddy Simulations of Stratified Atmospheric Boundary Layers.”

Organize.py and workflow improvements

Executing organize.py results in

425: Too Early: Seems like snakemake was never executed:

Either this should not be an error, or unlock should be fixed

Extend statistics toolbox

In order to calculate

• Wall shear stress, and
• SGS stress

The default statistics toolbox needs to be extended.

Frequency

1. Sampling frequency ~ inner units: u+, y+. Say every 10 to 20 time steps.
2. Averaging duration ~ eddy turn over time (1 day / 1 inertial oscillation)?

On-the-fly

Calculate <\tau_ij>, \tau_ij^2, <\tau_ij S_ij> .... (cf. SIMPSON documentation)

• Compute wall shear stress, u_rms, SGS stress
• Coordinate transformation (x,y,z) -> (x, z, y), informing the code which are the coordinates (wall normal, tangential)

Post-processing mode

• Interpolate to a regular grid
• Average in streamwise direction, to convert 2D statistics to 1D statistics

Considerations

1. Problem of eddy viscosity models, results in anisotropic stress
2. Needs to account for varying viscosity?
3. Ekman turning?

Plot horizontal energy spectra

Upload simplified package

For every push and for every tag.

Use https://github.com/actions/upload-artifact

Case log: Maronga et al. (2019)

• First try: 71eafc9, 9cb93f0
• Finer mesh resolution: 9cb93f0
• Filter weight 4, Filter cutoff 0.6: fab2528