This is an old revision of the document!
Table of Contents
Tornado Box
This page provides the Fortran program used for my 2009 QJ article. I also still offer a more ancient page that provides some results and animation, but no Fortran program.
The coding style evolved from a program passed down to Prof. Fiedler by his dissertation adviser, in the late 1970s. Notable ancient features are the Fortran 77 style, for the .f files, and the use of deprecated common blocks. Incidentally, the name of the array used to store all the parameters also has ancient roots: ap comes from the dissertation adviser.
Remarkably, the program still compiles using the Linux Fortran Compiler: ifx.
After you install the Fortran compiler, you need to execute this statement:
source /home/bfiedler/intel/oneapi/setvars.sh
and you should then see these lines in your terminal:
:: initializing oneAPI environment ... -bash: BASH_VERSION = 5.0.17(1)-release args: Using "$@" for setvars.sh arguments: :: compiler -- latest :: debugger -- latest :: mpi -- latest :: umf -- latest :: oneAPI environment initialized ::
Compilation should now be simple. In my shell program compileit you will see these
statements:
# assumes the following line already was executed: # source /home/bfiedler/intel/oneapi/setvars.sh ifx -o gridmaker gridmaker.f rm box.size ln -s box91.size box.size ifx -o boxn91 boxn.f boxs.f ifx -o dubl91 dubl.f ifx -o cubit cubit.f where.f90 # compile higher resolution model: rm box.size ln -s box181.size box.size ifx -o boxn181 boxn.f boxs.f
If you execute the compilieit.sh, in about 10 seconds you will have all the compiled Fortran programs that you will need. Then execute gridmaker. That will produce grid.dat and grid2d.dat. The latter was used in a companion axisymmetric model, that is not offered here. The file grid.dat is
a text file that contain the x,y,z coordinates of the stretched 3-D grid, of $91 \times 91 \times 46$ grid points. Grid sizes are of $181 \times 181 \times 91$ are also employed, but those grids are interpolated from the $91 \times 91 \times 46$, when the resolution is doubled in simulations to focus on events of interest.
test an animation from 2025:
And here is the domain, grid and buoyancy in the model, and also low-pressure isosurfaces of a vortex in the center:
Ancient, to be revised:
Assumes:
- You have the Intel Fortran compiler installed in
~/intel, and that you
untar the boxprogram directory in your home directory, meaning ~/boxprogram.
- You have the Anaconda distribution of python. At the time of this writing,
mayavi(for 3D renderings) only runs in python 2.7, and not python 3.0. Here is what I see
[L boxprogram]$ python Python 2.7.6 |Anaconda 2.0.0 (64-bit)| (default, May 27 2014, 14:50:58)
- you have
box.sizeand the header withingridmaker.sizeconfigured for the coarsest resolution:
PARAMETER(In=91,Jn=91,kn=46,
& im=In-1,Jm=Jn-1,km=kn-1)
$ cd ~ $ gunzip boxprogram.tar.gz $ tar xvf boxprogram $ cd boxprogram $ chmod u+x sourceintel $ ./sourceintel $ ifort -o gridmaker gridmaker.f $ ./gridmaker
Two files will be produced:
grid2d.datis for comparing with the companion axisymmetric simulations and is not needed.grid.datis what is needed for 3D simulations.
$ chmod u+x compileit $ ./compileit
That will produce two fortran programs:
boxn91The prognostic CFD programcubitInterpolates the fields within -.25<x<.25, -.25<y<.25 and 0<z<.5 in the output files fromboxn91, the.savfiles, to a regular Cartesian grid with equal spacing in all directions, as required bymayavi
Now we run the model boxn91.
$ cd ~ $ mkdir out_test $ ln -s out_test boxout $ cd boxprogram
Cases are identified with prefixes of four characters, which is short for historical reasons: it fits within four bytes. For example, aaaa.input configures case {aaaa} with the viscosity (inverse Reynolds number) to be rei=.0004 and the Coriolis parameter f to be swirl=.05. Somewhere within the simulation you should be able to get a snapshot like this, captured with the deprecated vis5d software:
running the model
The frequency at which the model dumps all the fields (and time derivatives of fields, so the model can in principle be exactly restarted) are within filetimes.txt.
$ ln -s aaaa.input box.input $ ./boxn91
This is what you enter at the prompts. It will take the model out to the dump of the dump of ~/boxout/aaaad010000.sav.
[L boxprogram]$ ./boxn91 continue from previous case? 0 enter four characters: 'case' 'aaaa' Current parameters are: &F REI = 3.9999999E-04, REIT = 1.0000000E-03, SWIRL = 5.0000001E-02, DTMAX = 0.1000000 , ALIM = 0.2000000 , DIFF = 0.1000000 , PIT = 3.000000 , RLX = 1.500000 , GAMMA = 0.3000000 , SLIP = 2.000000 , ORDER = 5.000000 , EX2 = 0.0000000E+00, EX3 = 0.0000000E+00, EX4 = 0.0000000E+00 / Do you want to change any parameters? enter yes, and store=2, yes, no store=1, no=0 0 dtdif= 2.0603456E-02 time now = 0.0000000E+00 enter tstop,tdiag 10,.01
After about ten minutes of wall clock, the model should make it out to time=10 (dimensionless).
4844. 9.980 0.002060 0.67087 0.67087 -0.04251
4849. 9.990 0.002060 0.67100 0.67100 -0.04268
4854. 10.001 0.002060 0.67112 0.67112 -0.04284
FINISHED!
nloc= 4854 elapsed/nloc 0.0000000E+00 millisecs
time= 10.00070 cput= 0.0000000E+00 steps= 4854.000
The numbers you see are also within aaaa.hist, which will be appended to after restart. Here is what the numbers mean:
4854. 10.001 0.002060 0.67112 0.67112 -0.04284
step # time dt max speed max w min pressure
The numbers are diagnosed for the central part of the domain, and z<.5, roughly the domain portion seen in the above snapshot.
Next we show how to restart the model.
$ rm restart.sav $ ln -s ../boxout/aaaad010000.sav restart.sav $ ./boxn91
Here is way it looks for me. Note there is a dump of the vertical profile of viscosity. Apparently I was having some confusion/problem about this at some time, so that is why the dump shows up. Note a the number 20.00000 is spit out, which is the time for the next file dump, as acquired from filetimes.txt.
[L boxprogram]$ ./boxn91
continue from previous case?
1
1 0.0000000E+00 3.9999999E-04 3.9999999E-04
2 4.4620000E-03 3.9999999E-04 3.9999999E-04
3 9.0290001E-03 3.9999999E-04 3.9999999E-04
4 1.3805000E-02 3.9999999E-04 3.9999999E-04
5 1.8893000E-02 3.9999999E-04 3.9999999E-04
6 2.4390001E-02 3.9999999E-04 3.9999999E-04
7 3.0393001E-02 3.9999999E-04 3.9999999E-04
8 3.6991000E-02 3.9999999E-04 3.9999999E-04
9 4.4270001E-02 3.9999999E-04 3.9999999E-04
10 5.2308001E-02 3.9999999E-04 3.9999999E-04
11 6.1175998E-02 3.9999999E-04 3.9999999E-04
12 7.0941001E-02 3.9999999E-04 3.9999999E-04
13 8.1660002E-02 3.9999999E-04 3.9999999E-04
14 9.3382001E-02 3.9999999E-04 3.9999999E-04
15 0.1061500 3.9999999E-04 3.9999999E-04
16 0.1200000 3.9999999E-04 3.9999999E-04
17 0.1349580 3.9999999E-04 3.9999999E-04
18 0.1510460 3.9999999E-04 3.9999999E-04
19 0.1682760 3.9999999E-04 3.9999999E-04
20 0.1866550 3.9999999E-04 3.9999999E-04
21 0.2061860 3.9999999E-04 3.9999999E-04
22 0.2268610 3.9999999E-04 3.9999999E-04
23 0.2486730 3.9999999E-04 3.9999999E-04
24 0.2716050 3.9999999E-04 3.9999999E-04
25 0.2956400 3.9999999E-04 3.9999999E-04
26 0.3207550 3.9999999E-04 3.9999999E-04
27 0.3469230 3.9999999E-04 3.9999999E-04
28 0.3741180 3.9999999E-04 3.9999999E-04
29 0.4023070 3.9999999E-04 3.9999999E-04
30 0.4314580 3.9999999E-04 3.9999999E-04
31 0.4615390 3.9999999E-04 3.9999999E-04
32 0.4925120 3.9999999E-04 3.9999999E-04
33 0.5243420 3.9999999E-04 3.9999999E-04
34 0.5569940 3.9999999E-04 3.9999999E-04
35 0.5904310 3.9999999E-04 3.9999999E-04
36 0.6246150 3.9999999E-04 3.9999999E-04
37 0.6595120 3.9999999E-04 3.9999999E-04
38 0.6950850 3.9999999E-04 3.9999999E-04
39 0.7313000 4.6260000E-04 4.6260000E-04
40 0.7681220 5.3624407E-04 5.3624407E-04
41 0.8055170 6.1103405E-04 6.1103405E-04
42 0.8434540 6.8690802E-04 6.8690802E-04
43 0.8819000 7.6380005E-04 7.6380005E-04
44 0.9208260 8.4165210E-04 8.4165210E-04
45 0.9602020 9.2040398E-04 9.2040398E-04
46 1.000000 1.0000000E-03 1.0000000E-03
continuing with aaaa
20.00000
Current parameters are:
&F
REI = 3.9999999E-04,
REIT = 1.0000000E-03,
SWIRL = 5.0000001E-02,
DTMAX = 0.1000000 ,
ALIM = 0.2000000 ,
DIFF = 0.1000000 ,
PIT = 3.000000 ,
RLX = 1.500000 ,
GAMMA = 0.3000000 ,
SLIP = 2.000000 ,
ORDER = 5.000000 ,
EX2 = 0.0000000E+00,
EX3 = 0.0000000E+00,
EX4 = 0.0000000E+00
/
Do you want to change any parameters?
enter yes, and store=2, yes, no store=1, no=0
0
dtdif= 2.0603456E-02
time now = 10.00070 enter tstop,tdiag
200,.01
You will need more than 190 minutes for the model to complete. As the tornado spins up, the model reduces the time step to satisfy the CFL condition.
plotting some results
After the model finishes, you should have time series in aaaa.hist that can be plotted.
$ python plothist.py aaaa
You should see something like this, traces of the maximum values of speed, w and q:
Note q in the plot is a way of plotting pressure as a speed. It is the speed that Bernuolli's equation predicts could be accelerated by the low pressure $q=\sqrt{2|p|}$.
The absolute value just keeps the code from blowing up; the negative values of p are producing the extreme values of q that are plotted.
pickle files
The python plotting programs for the fields within the .sav files require first converting the .sav files to a pickle file. This one-time conversion makes it very easy to read the data into a variety of python plotting programs. Here we do an example with the final output at time=200.
$ python sav2pkl.py ~/boxout/aaaad020000.sav
To see the grid at the surface and in a vertical cross-section through the center:
$ python gridplot.py ~/boxout/aaaad020000.pkl $ display grid_xy.svg $ display grid_xz.svg
To a see an ugly plot of the vertical velocity at z=.06, together with horizontal wind vector plotted on 1/16th of the grid points, do:
$ python plotpkl.py ~/boxout/aaaad020000.pkl $ eog w20000.png
Plots of the fields will be more effective after we focus on the region where the vortex is. The domain is -2<x<2 -2<y<2 and 0<z<1. We make a new .sav file that stores cuts out -.25<x<.25 -.25<y<.25 and 0<z<.5 and interpolates the fields to a Cartesian, isotropic grid.
Currently, cubit is a bit clumsy to run. Here is how we do it:
$ cd ../boxout $ ~/boxprogram/cubit
You need to enter some text at the query:
[L out_aaaa]$ ~/boxprogram/cubit enter 'case' ,num 'aaaa', 200000
Here is the tail end of what I see dumped to the monitor:
interpolated grid size= 3.9062500E-03 minimum model grid size= 1.3805000E-02 aaaad200000.sav aaaaa200000.sav
Note the fifth character in the original .sav file is d, the altered (interpolated) .sav file has an a. There is no particularly good reason for this designation.
$ cd ~/boxprogram $ python sav2pkl.py ~/boxout/aaaaa200000.sav $ python plotpkl.py ~/boxout/aaaaa200000.pkl $ eog w20000.png
Here is what I see:
3d renderings with mayavi
python isovectpkl.py ~/boxout/aaaaa200000.pkl
A GUI should pop open, showing the p=-1 and p=-.5 isosurfaces. Note the upper left icon allows some features of the rendering to be switched on or off. Dexerity and creativity with a two-button mouse is also essential to change the viewing of the rendering.
doubling the resolution
Note that it is possible to identify a time period of interest and restart the model from a .sav file, with more dump times specified within filetimes.txt.
It is also possible to double the resolution before restarting.
$ cd ~/boxout $ ../boxprogram/dubl181
Here is what I did:
[L out_aaaa]$ ../boxprogram/dubl181 enter 'case' ,num, 'newcase' 'aaaa',200000,'bbbb'
This makes a new file bbbbE200000.sav.. By the way, the E does not mean much anything special … maybe Extended? In order to start from this file, you need to edit box.size to be
c PARAMETER(In=361,Jn=361,kn=181,
c PARAMETER(In=91,Jn=91,kn=46,
PARAMETER(In=181,Jn=181,kn=91,
& im=In-1,Jm=Jn-1,km=kn-1)
and then compile again
$ ifort -o boxn181 boxn.f boxs.f
Proceed with your computations using boxn181
left undone
What isn't shown in this tutorial is how to make animations, such as seen here: http://mensch.org/suction/
box (last edited 2014-06-17 10:24:43 by BrianFiedler)