# Tutorials¶

To run HERCULES, you need the main program **hercules.exe** (which should be generated after the installation) and the input file **parameters.input**. See the detailed explanation of the input parameters at the end of this page. In addition, you need to put the initial turbulent flow fields in the **results** folder. Otherwise, it will take a very long time to generate turbulence from a laminar flow field.

When you run **hercules.exe**, the code will try to read initial fields named **init_u.dat, init_v.dat**, etc. from the **results** folder. The best way to get the initial flow fields is from your previous simulations. Note that HERCULES will save intermediate flow fields named **bak_u.dat, bak_v.dat**, etc. for re-start runs. You can rename them to **init_u.dat, init_v.dat**, etc.

The above trick will not work if you are running HERCULES for the first time. Therefore, I prepare some coarse-resolution initial turbulent flow fields (e.g., **coarse_u_32_cubic.dat**) for your first run. These coarse flow fields come with HERCULES, they are in the **misc/Turb_Gen** folder. To interpolate these coarse-resolution flow fields to any resolution you need, you can use the Matlab code **Interpolate_Initial_Fields.m** provided in the **misc/Turb_Gen** folder. After you run the Matlab code, some initial fields named **init_u.dat, init_v.dat**, etc. will be generated. You need to move these files to the **results** folder.

**NOTE**: if you don’t have an access to Matlab, use Octave (https://www.gnu.org/software/octave) instead. Octave is a free alternative to Matlab and it can be installed on Ubuntu by running:

```
sudo apt-get install octave
```

To run the Matlab script using Octave, do:

```
octave Interpolate_Initial_Fields.m
```

The followings are the specific steps you need to follow to run HERCULES for a few DNS benchmarks.

**Plane Closed-channel Flows at Re_tau=180**This is a classical DNS benchmark case for turbulent channel flows, please refer to Moser, Kim, and Mansour (1999), “Direct numerical simulation of turbulent channel flow up to Re_tau=590” in Physics of Fluids. (hereinafter referred to as MKM99).

Since the default setup of HERCULES is for MKM99, you don’t need to change

**parameters.input**. You only need to prepare some initial turbulent fields. Here are the steps you need to follow.Go to the

**misc/Turb_Gen**folder, and run the Matlab script**Interpolate_Initial_Fields.m**. This script will generate some files named**init_u.dat, init_v.dat**, etc.Move the generated files

**init_u.dat, init_v.dat**, etc. to the**results**folder.In the main folder, start the simulation using 2 CPU cores by running:

mpirun -np 2 hercules.exe

The simulation results (including the instantaneous, mean, and re-start files) will be stored in the

**results**folder. These files are in binary format. To plot the results (e.g., the mean profiles), you can use the Matlab script**Plot_Mean_Profiles.m**in the**/misc/Post_Processing**folder.By default, the simulation will be run for 20 non-dimensional times. You probably need a longer time for a fully developed turbulent field. To do this, after the above simulation is done, you can rename the re-start files

**bak_u.dat**to**init_u.dat**(same for other variables). Then, you can re-run the simulation with the new initial fields:mpirun -np 2 hercules.exe

If you want to use more CPU cores, you need to change p_row and p_col in

**parameters.input**accordingly. For example, you want to run HERCULES using 8 cores, you need to change p_row=4 and p_col=2 in**parameters.input**, and run:mpirun -np 8 hercules.exe

So far, you should be able to run HERCULES for the MKM99 case. Good luck!

**Stably Stratified Closed-channel Flows at Re_tau=180 and Ri_tau=18**This is a benchmark case for stratified turbulent channel flows, please refer to Garcia-Villalba and del Alamo (2011), “Turbulence modification by stable stratification in channel flow” in Physics of Fluids (hereinafter referred to as GD11).

To run HERCULES for GD11, you need to change some parameters in

**parameters.input**. In addition, you need to prepare some initial turbulent fields. Here are the steps you need to follow.Copy

**tutorials/parameters.input.GD11**to the main folder and rename it to**parameters.input**; replace the old**parameters.input**.Go to the

**misc/Turb_Gen**folder. In the Matlab script**Interpolate_Initial_Fields.m**, change nx2=288, ny2=288, nz2=128. These are the grid numbers we need for the GD11 case. Run the script, and it will generate some files named**init_u.dat, init_v.dat**, etc.Move the generated files

**init_u.dat, init_v.dat**, etc. to the**results**folder.In the main folder, start the simulation using 2 CPU cores by running:

mpirun -np 2 hercules.exe

You can follow the steps 4 and 5 in the MKM99 case for plotting the output and doing the re-start runs.

**Neutrally Stratified Ekman Layer Flows at Re_g=400**This is a benchmark case for Ekman layer flows, please refer to Coleman, Ferziger, and Spalart (1990), “A numerical study of the turbulent Ekman layer” in Journal of Fluid Mechanics. (hereinafter referred to as CFS90). Here are the steps you need to follow.

Copy

**tutorials/parameters.input.CFS90**to the main folder and rename it to**parameters.input**; replace the old**parameters.input**.Go to the

**misc/Turb_Gen**folder. In the Matlab script**nterpolate_Initial_Fields.m**, change nx2=288, ny2=288, nz2=128. In addition, change Uscaling=0.06. Note that we need to scale the initial field for the CFS90 case. Run the script, and it will generate some files named**init_u.dat, init_v.dat**, etc.Move the generated files

**init_u.dat, init_v.dat**, etc. to the**results**folder.In the main folder, start the simulation using 2 CPU cores by running:

mpirun -np 2 hercules.exe

Again, the CFS90 case needs a longer time to converge, you can follow the steps 4 and 5 in the MKM99 case for plotting the output and doing the re-start runs.

Here is a sample

**parameters.input**file for the MKM99 case. The meaning of each parameter is explained as follows:&domain p_row = 2 p_col = 1 lx = 12.56637061 ly = 4.188790205 lz = 2.0 nx = 256 ny = 192 nz = 128 zstretch = 1.1 ochannel = 0 / &modeling dp_opt = 1 cds = 2 issk = 1 dts = 2 isnoise = 0 noise_mag = 0.1 isdamp = 0 nzdamp = 15 cdamp = 0.01 / &constants nu = 0.005555556 alpha = 0.007936508 got = 0.0 fc = 0.0 ug = 0.0 vg = 0.0 dt = 0.001 imax = 20000 u_mrf = 0.0 isscalar = 0 t_ref = 0 tbot = -0.5 ttop = 0.5 is_ri_var = 0 ri_str = 0.0 ri_end = 0.0 / &io restart = 1 istmsr = 1 ibackup = 2000 iinstfl = 2000 imeanfl = 2000 isxy2d = 0 xy2d_id = 10,32,64 isxz2d = 0 xz2d_id = 32,64 isyz2d = 0 yz2d_id = 32,64 intv_2d = 1000 /

**p_row, p_col**- These two parameters control the domain decomposition in the x and y directions for parallel simulation. The grid numbers nx and ny should be divisible by p_row, and the grid numbers ny and nz should be divisible by p_col. Make sure to use p_row*p_col cores for parallel simulation. To ensure high parallel efficiency, p_row shouldn’t be too close to nx and ny, and p_col shouldn’t be too close to ny and nz. In addition, p_row and p_col should be as close as possible.
**lx, ly, lz**- Domain sizes in the x, y, and z directions.
**nx, ny, nz**- Number of cells in the x, y, and z directions.
**zstretch**- This parameter controls the grid refinement near the wall. The larger zstretch the finner grids near the wall. The default value for zstretch is 1.1.
**ochannel**- This parameter controls the channel configuration. 0: closed channel flows (non-slip boundary condition for the top and bottom walls). 1: open channel flows (non-slip and slip boundary conditions for the bottom and top walls, respectively).
**dp_opt**- This parameter controls the options to drive the flow. 1: constant friction. Re and Ri are friction Reynolds number and Richardson number in this case. 2: constant bulk velocity. Re and Ri are bulk Reynolds number and Richardson number in this case. 3: constant geostrophic wind for Ekman layer. Re and Ri are geostrophic Reynolds number and and bulk Richardson number in this case.
**cds**- Differential scheme for spatial derivatives in the horizontal directions. 1: 2nd order central difference scheme. 2: 4th order central difference scheme. 3: Spectral method.
**issk**- If the nonlinear terms are casted in skew-symmetric form (only for cds=3). For the spectral method, the skew-symmetric form has a better numerical stability (compared with the divergence form) and is strongly recommended. 0: no. 1: yes.
**dts**- Time advancement scheme for the convective terms (for the diffusion term, the Crank-Nicolson scheme is used.) 1: 2nd order Adam-Bashforth. 2: 3rd order Runge-Kutta.
**isnoise, noise_mag**- Add random noise for the initial fields (to generate turbulence). isnoise=1: add random noise to the initial fields. noise_mag: the magnitude (%) of the random noise if isnoise=1.
**isdamp, nzdamp, cdamp**- Rayleigh damping parameters. isdamp=1: apply Rayleigh damping to the domain top. nzdamp: how many points for damping if isdamp=1. cdamp: damping coefficient if isdamp=1.
**nu**- Kinematic viscosity. For the constant friction simulations (dp_opt=1), set nu=1/Re_tau. Similarly, for the constant mass flow rate simulation, set nu=1/Re_b. For the Ekman layer case, set nu=1/Re_g.
**alpha**- Thermal diffusivity. Usually, we set alpha=nu/Pr, where Pr=0.7.
**got**- g/T. Here g is the gravitational acceleration and T is temperature. Set this parameter to non-zero for stratified flow simulation. For the constant friction simulations (dp_opt=1), set got=Ri_tau. Similarly, for the constant mass flow rate simulation, set got=Ri_b. For the Ekman layer case, set got=Ri_g.
**fc**- Coriolis parameter. This parameter is for Ekman layer simulation. Set it to nu*2.0.

**ug, vg**
Geostrophic wind speeds. Set ug to 1.0 for Ekman layer simulation.

**dt**- Time step. dt is restricted by the maximal CFL number.
**imax**- How many steps to run.
**u_mrf**- Moving speed of the reference frame. The default value is 0. Set it to be about half of the bulk velocity can reduce the maximal CFL number.
**isscalar**- Whether to include scalar (temperature) for the simulation.
**tbot and ttop**- Bottom wall and top wall temperature. For stratified closed channel flow case, set tbot=-0.5 and ttop=0.5. For stratified open-channel case, set tbot=-1.0 and ttop=0.0.
**is_ri_var and ri_str and ri_end**- These parameters define the time-varying Richardson number value during the simulation. is_ri_var=1: Richardson number will vary during the simulation. ri_str: initial ri value if is_ri_var=1. ri_end: final ri value if is_ri_var=1.
**restart**- 1: read initial field from “results/init_u.dat”, “results/init_v.dat”, etc. 0: use a default initial field (u=0,v=0,w=0,t=0,p=0).
**istmsr**- istmsr=1: output time-series of u, v, w, p, t collected from a vertical line in the center of the domain.
**ibackup, iinstfl, imeanfl**- These parameters control the output frequency of the backup, instantaneous, and mean fields. ibackup: output frequency of the backup fields. iinstfl: output frequency of the instantaneous fields. imeanfl: output frequency of the mean fields.
**isxy2d, xy2d_id**- isxy2d=1: output 2D slices of u, v, w, t at specific x-y planes. xy2d_id: the k index for the 2D slice, you can give values for multiple levels.
**isxz2d,xz2d_id**- isxz2d=1: output 2D slices of u, v, w, t at specific x-z planes. xz2d_id: the j index for the 2D slice, you can give values for multiple levels.
**isyz2d, yz2d_id**- isyz2d=1: output 2D slices of u, v, w, t at specific y-z planes. yz2d_id: the i index for the 2D slice, you can give values for multiple levels.
**intv_2d**- The output frequency of the 2D slices output (steps).