Thoughts

Grab the code: tornadobox at gitbub.

This page describes 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 downloading, execute the shell program:

$ ./intel-fortran-compiler-2025.0.4.21.sh

Then do:

$ source ~/intel/oneapi/setvars.sh

You should then see in your terminal a few lines:

:: 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 ::

Then, as per instructions, enter into hello.f90:

print *, "Hello, world!"
end

Then compile:

$ ifx -o hello hello.f90

That might not have worked. When I first tried ifx, I maybe saw this:

ld: cannot find -lm: No such file or directory
ld: cannot find -lpthread: No such file or directory
ld: cannot find -lc: No such file or directory

Bummer…

Then I think I did this (as suggested by chatGPT ?):

sudo apt install gfortran

What gfortran has to do with ifx, I don't know. Miraculously, hello.f90 then compiled, and ./a.out executed the test program.

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.

The program cubit was used in the past as a needed precursor to interpolating to the Cartesian grids that were required by the deprecated Vis5d, which had been used for the renderings. Here we used python and the module pyvista for the renderings.

Here is the domain, grid and buoyancy in the model, and also low-pressure isosurfaces of a vortex in the center at time t=15, for $\mu=.0001$ and $f=.1$:

The gaussian fixed buoyancy field has ellipsoid isosurfaces, three of which are shown. The buoyancy could sustain a hydrostatic pressure deficit of 0.5 at the surface. Alternatively the buoyancy could accelerate a parcel to velocity of 1. The pressure isosurfaces are -0.5 and -1.0. The vertical velocity in the vortex core is $w=1.43$.

And here we continue with the simulation from $t=15$ to $t=30$, after doubling the resolution, and zooming in, with a video:

Model runs are identified with a four-character case prefix. Here we suggest the first letter is a, b or c for low, medium or high value of viscosity $\mu$. Then $f$ and two digits to denote swirl, which is essentially a Coriolis parameter. If the case is high resolution, we use a capital $F$. We offer bf10 and bF10 in the examples here.

Here is a plot of the time history of extreme values of speed $s$, vertical velocity $w$, and $q$. Note $q=\sqrt{2|p|}$ in the plot is a way of plotting pressure as a speed. $q$ can be interpreted as the speed that Bernoulli's equation predicts could be accelerated by the low pressure .

The FORTRAN programs do not accept command line arguments. They may read in files and, when started, need a few strings to be input in response to queries. When first started, filetimes.txt and box.input must be in the current working directory. boxn91 and boxn181 will then output the fields from the simulation at the times specified in filetimes.txt. These .sav files are output in a directory ../boxout relative to the current working directory. It is highly recommended that you create a directory with the 4-character case name, and then create a symbolic link boxout to that directory. Also, box.input is recommended to be a symbolic link. An input file bf10.input is included with the distro. If you want to explore the parameter space, you can copy bf10.input to a new named .input file, and edit it. For my research, I only change the inverse Reynolds number REI, or $\mu$, and SWIRL, which is $f$.

Let us begin.

$ cd tornadobox
$ mkdir ../bf10
$ ln -s bf10 ../boxout
$ ./gridmaker
$ ./makefiletimes.py 0 202 1
$ ./makefiletimes.py 0 202 1 > filetimes.txt
$ ./boxn91
 continue from previous case?
0
  enter four characters: 'case'
'bf10'
  Current parameters are:
 &F
 REI     =  9.9999997E-05,
 REIT    =  1.0000000E-03,
 SWIRL   =  0.1000000    ,
 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
 time now =  0.0000000E+00  enter tstop,tdiag
2,.01

Note: I always answer 0 to the question Do you… In less than a minute the model will say FINISHED!

Now let's restart the model. But first, put this line at then end of your .bashrc:

export PATH="/home/bfiedler/tornadobox:/home/bfiedler/bin:$PATH"

And, if you haven't done so, you may want this also:

source /home/bfiedler/intel/oneapi/setvars.sh

Then open a new window, and tornadobox should now be in your path.

$ cd ~/bf10
$ ls
$ cat bf10.hist
$ ln -s bf10d002000.sav restart.sav
$ makefiletimes.py 0 202 1 > filetimes.txt
$ boxn91
 continue from previous case?
1
  continuing with bf10
  time now:    2.002613    
 will write .sav to disk at time:   3.000000    
  Current parameters are:
 &F
 REI     =  9.9999997E-05,
 REIT    =  1.0000000E-03,
 SWIRL   =  0.1000000    ,
 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
 time now =   2.002613      enter tstop,tdiag
5,.01

Note: the answer 1 to continue from previous case? </code>

Let's restart the model again, this time so that we can logout from the interactive session. Use a text editor to make a file name answers, and put in the 3 lines:

1
0
30,.01

Then

$ remove restart.sav
$ ln -s bf10d005000.sav restart.sav
$ nohup boxn91 < answers > outstuff.txt &
$ tail -f outstuff.txt
$ ^C
$ logout

Let's have a break. When you come back check to see FINISHED in outstuff.txt. On one of my ancient computers with Intel(R) Core(TM) i5-6500 CPU @ 3.20GHz, that required 45 minutes of wall clock.

double the resolution, and repeat the simulation from $t=15$ to $t=30$.

$ dubl91 bf10d015000.sav
  enter 'case' ,num, 'newcase' 
'bf10',015000,'bF10'

$ ln -s bF10E015000.sav restart.sav $ makefiletimes.py 0 202 .05 > filetimes.txt

continue from previous case? 1

continuing with bF10
time now:    15.00064    

will write .sav to disk at time: 15.05000

Current parameters are:

&F REI = 9.9999997E-05, REIT = 1.0000000E-03, SWIRL = 0.1000000 , 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 time now = 15.00064 enter tstop,tdiag 15.1,.01 </code> That should finish in a few minutes. To simulate out to $t=30$, I needed about 11 hours and 100G of disk space. Your previous answers file should still be appropriate.

$ rm restart.sav
$ ln -s bF10D015100.sav restart.sav
$ nohup boxn181 < answers >outstuff.txt &
$ tail -f outstuff.txt
$ ^C
$ logout

Check it in the morning.

Let's make the three images seen above. Again, assuming your tornadobox directory is in your path:

$ cd ~/boxout
$ hist2png.py  bf10.hist bF10.hist

Now let's plot the isosurfaces of the fields.

$ ls b*.sav

You should see you have a lot of .sav files. Now lets extract what we need into python pickle files. Such files are easily readable by our plotting programs. This may need about 20 minutes:

$ sav2pkl.py b*.sav

If your python does not have the pyvista module, you will need to install it. To make the times tidy and complete, let's rename a file:

$ mv bF10E015000.pkl bF10D015000.pkl
$ pkl2png.py -s -hg -vg   bF10D01500

Note that pkl2png.py will do the glob of bF10D01500*.pkl, in this case finding one .pkl file. The -s option should have shown the image do the screen, like the first one at the top of this page.

Next we zoom in to 25% of the camera distance, and plot at $t=30$:

$ pkl2png.py -s -hg -c .25 bF10D030000

And add a $w$ isosurface:

$ pkl2png.py -s -hg -w -c .25 bF10D030000

$ Now we go for it, and plot the frames for movie. No -s option:

$ mkdir pngs
$ pkl2png.py -hg -w -c .25 bF10D0 pngs

That required 10 minutes on my computer to produce 301 .png files.

$ cd pngs
$ eog bF10D015000.png

The eog image viewer allows for the right arrow key to rapidly step through the images.

$ cd ~/boxout
$ maketitle.py pngs/bF10D015000.png 
$ pip install moviepy
$ png2mp4.py -t title.png pngs/bF10D
$ mpv video.mp4

FINISHED!