Thoughts

User Tools

Site Tools

Trace:
tornadobox

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
tornadobox [2025/03/18 16:22] – admintornadobox [2025/04/13 07:58] (current) – admin
Line 1: Line 1:
-==== Tornado Box====+==== 3-D Fortran CFD code for tornado-like vortices ====
  
-This page provides the Fortran program used for +Grab the code: [[https://github.com/bfiedler/tornadobox | tornadobox at gitbub]]. 
 + 
 +This page describes the Fortran program used for 
 {{ :tornadobox:suction_vortices_and_spiral_breakdown_in_numerical.pdf | my 2009 QJ article}}. {{ :tornadobox:suction_vortices_and_spiral_breakdown_in_numerical.pdf | my 2009 QJ article}}.
 I also still offer [[https://12characters.net/suction/|a more ancient page]] I also still offer [[https://12characters.net/suction/|a more ancient page]]
 that provides some results and animation, but no Fortran program. 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. +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. Feel free to  
 +see the Fortran code at [[https://12characters.net/tornadoboxbrowser/ | tornadobox browser]]. 
  
 === FORTRAN === === FORTRAN ===
Line 75: Line 78:
  
 The program ''cubit'' was used in the past as a needed precursor to interpolating to the Cartesian grids that were required by the deprecated [[https://en.wikipedia.org/wiki/Vis5D | Vis5d]], which had been used for the renderings.  Here we used ''python'' and the module ''pyvista'' for the renderings.  The program ''cubit'' was used in the past as a needed precursor to interpolating to the Cartesian grids that were required by the deprecated [[https://en.wikipedia.org/wiki/Vis5D | Vis5d]], which had been used for the renderings.  Here we used ''python'' and the module ''pyvista'' for the renderings. 
 +
 +Make your own browser with [[https://github.com/bfiedler/f90tohtml|f90tohtml]]:
 +<code>
 +f90tohtml tornadobox.f2h
 +</code>
  
 === An example case === === An example case ===
Line 206: Line 214:
 $ ^C $ ^C
 $ logout $ logout
- 
-Let's have a break. 
- 
-=== ancient, to be revised === 
-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  
- 
- 
- 
-<code> 
-[L boxprogram]$ python 
-Python 2.7.6 |Anaconda 2.0.0 (64-bit)| (default, May 27 2014, 14:50:58)  
 </code> </code>
  
 +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$.
-  * you have ''%%box.size%%'' and the header within ''%%gridmaker.size%%'' configured for the coarsest resolution:  +
- +
  
 <code> <code>
-        PARAMETER(In=91,Jn=91,kn=46, +$ dubl91 bf10d015000.sav 
-     &            im=In-1,Jm=Jn-1,km=kn-1) +  enter 'case' ,num, 'newcase'  
-</code>+'bf10',015000,'bF10' 
 +$ ln -s bF10E015000.sav restart.sav 
 +$ makefiletimes.py 0 202 .05 > filetimes.txt
  
- +continue from previous case? 
- +1 
- +  continuing with bF10 
-<code> +  time now:    15.00064     
-$ cd ~ + will write .sav to disk at time:   15.05000    
-$ gunzip boxprogram.tar.gz +
-$ tar xvf boxprogram +
-$ cd boxprogram +
-$ chmod u+x sourceintel +
-$ ./sourceintel +
-$ ifort -o gridmaker gridmaker.f +
-$ ./gridmaker +
-</code> +
- +
- +
-Two files will be produced:  +
- +
-  * ''%%grid2d.dat%%'' is for comparing with the companion axisymmetric simulations and is not needed.  +
-  * ''%%grid.dat%%'' is what is needed for 3D simulations.  +
- +
- +
- +
-<code> +
-$ chmod u+x compileit +
-$ ./compileit +
-</code> +
- +
- +
-That will produce two fortran programs:  +
- +
-  * ''%%boxn91%%'' The prognostic CFD program  +
-  * ''%%cubit%%'' Interpolates the fields within -.25<x<.25, -.25<y<.25 and 0<z<.5 in the output files from ''%%boxn91%%'', the ''%%.sav%%'' files, to a regular Cartesian grid with equal spacing in all directions, as required by ''%%mayavi%%''  +
- +
-Now we run the model ''%%boxn91%%''.  +
- +
- +
- +
-<code> +
-$ cd ~ +
-$ mkdir out_test +
-$ ln -s out_test boxout +
-$ cd boxprogram +
-</code> +
- +
- +
-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:  +
- +
-{{:f0p05-reip0004.png |}} +
- +
- +
- +
-===== 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%%''.  +
- +
- +
- +
-<code> +
-$ ln -s aaaa.input box.input +
-$ ./boxn91 +
-</code> +
- +
- +
-This is what you enter at the prompts. It will take the model out to the dump of the dump of ''%%~/boxout/aaaad010000.sav%%''.  +
- +
-<code> +
-[L boxprogram]$ ./boxn91 +
- continue from previous case? +
-0 +
-  enter four characters: 'case' +
-'aaaa'+
   Current parameters are:   Current parameters are:
  &F  &F
- REI     =  3.9999999E-04,+ REI     =  9.9999997E-05,
  REIT    =  1.0000000E-03,  REIT    =  1.0000000E-03,
- SWIRL   =  5.0000001E-02,+ SWIRL   =  0.1000000    ,
  DTMAX   =  0.1000000    ,  DTMAX   =  0.1000000    ,
  ALIM    =  0.2000000    ,  ALIM    =  0.2000000    ,
Line 330: Line 254:
   enter yes, and store=2, yes, no store=1, no=0   enter yes, and store=2, yes, no store=1, no=0
 0 0
- dtdif=  2.0603456E-02 + time now =   15.00064      enter tstop,tdiag 
- time now =  0.0000000E+00  enter tstop,tdiag +15.1,.01
-10,.01+
 </code> </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.
-After about ten minutes of wall clock, the model should make it out to ''%%time=10%%'' (dimensionless).  +
- +
- +
- +
-<code> +
-    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    +
-</code> +
- +
- +
-The numbers you see are also within ''%%aaaa.hist%%'', which will be appended to after restart. Here is what the numbers mean:  +
- +
-<code> +
-    4854.   10.001 0.002060   0.67112   0.67112  -0.04284 +
-    step #   time    dt       max speed  max w    min pressure +
-</code> +
- +
- +
-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.  +
 <code> <code>
 $ rm restart.sav $ rm restart.sav
-$ ln -s ../boxout/aaaad010000.sav restart.sav +$ ln -s bF10D015100.sav restart.sav 
-$ ./boxn91+$ nohup boxn181 < answers >outstuff.txt & 
 +$ tail -f outstuff.txt 
 +$ ^C 
 +$ logout
 </code> </code>
  
 +Check it in the morning.
  
-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%%''. +=== plotting results ===
  
 +Let's make the three images seen above.
 +Again, assuming your ''tornadobox'' directory is in your path:
 <code> <code>
-[L boxprogram]$ ./boxn91 +$ cd ~/boxout 
- continue from previous case? +$ hist2png.py  bf10.hist bF10.hist
-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+
 </code> </code>
  
- +Now let's plot the isosurfaces of the fields.
-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.  +
 <code> <code>
-$ python plothist.py aaaa+$ ls b*.sav
 </code> </code>
- +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:
-You should see something like this, traces of the maximum values of speed, w and q:  +
- +
- +
-{{:aaaa.png}}  +
- +
-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.  +
- +
- +
 <code> <code>
-$ python sav2pkl.py ~/boxout/aaaad020000.sav+$ sav2pkl.py b*.sav
 </code> </code>
- +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:
-To see the grid at the surface and in a vertical cross-section through the center:  +
 <code> <code>
-$ python gridplot.py ~/boxout/aaaad020000.pkl +$ mv bF10E015000.pkl bF10D015000.pkl 
-$ display grid_xy.svg +$ pkl2png.py -s -hg -vg   bF10D01500
-$ display grid_xz.svg+
 </code> </code>
 +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$:
-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:  +
- +
- +
 <code> <code>
-$ python plotpkl.py ~/boxout/aaaad020000.pkl +$ pkl2png.py -s -hg -c .25 bF10D030000
-$ eog w20000.png +
 </code> </code>
- +And add a $w$ isosurface:
- +
-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:  +
- +
- +
 <code> <code>
-$ cd ../boxout +$ pkl2png.py -s -hg -w -c .25 bF10D030000
-$ ~/boxprogram/cubit+
 </code> </code>
- +$ Now we go for it, and plot the frames for movie.  No ''-s'' option:
- +
- +
-You need to enter some text at the query:  +
- +
- +
 <code> <code>
-[L out_aaaa]$ ~/boxprogram/cubit +$ mkdir pngs 
-  enter 'case' ,num +$ pkl2png.py -hg -w -c .25 bF10D0 pngs
-'aaaa', 200000+
 </code> </code>
- +That required 10 minutes on my computer to produce 301 ''.png'' files.
- +
-Here is the tail end of what I see dumped to the monitor:  +
 <code> <code>
- interpolated grid size=  3.9062500E-03 +$ cd pngs 
- minimum model grid size=  1.3805000E-02 +$ eog bF10D015000.png
- aaaad200000.sav  aaaaa200000.sav+
 </code> </code>
- +The ''eog'' image viewer allows for the right arrow key to rapidly step through the images.
- +
-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.  +
- +
- +
- +
-<code> +
-$ cd ~/boxprogram +
-$ python sav2pkl.py ~/boxout/aaaaa200000.sav +
-$ python plotpkl.py ~/boxout/aaaaa200000.pkl +
-$ eog w20000.png  +
-</code> +
- +
- +
-Here is what I see:  +
- +
-{{:w20000.png}}  +
- +
-==== 3d renderings with mayavi ==== +
- +
- +
- +
- +
-<code> +
-python isovectpkl.py ~/boxout/aaaaa200000.pkl  +
-</code> +
- +
- +
-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.  +
- +
-{{:mayavi200.png}}  +
- +
-==== 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.  +
- +
- +
 <code> <code>
 $ cd ~/boxout $ cd ~/boxout
-$ ../boxprogram/dubl181+$ maketitle.py pngs/bF10D015000.png  
 +$ pip install moviepy 
 +$ png2mp4.py -t title.png pngs/bF10D 
 +$ mpv video.mp4
 </code> </code>
- +**FINISHED!**
- +
-Here is what I did:  +
- +
- +
- +
-<code> +
-[L out_aaaa]$ ../boxprogram/dubl181 +
-  enter 'case' ,num, 'newcase'  +
-'aaaa',200000,'bbbb' +
-</code> +
- +
- +
-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  +
- +
-<code> +
-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) +
-</code> +
- +
- +
-and then compile again  +
- +
- +
- +
-<code> +
-$ ifort -o boxn181 boxn.f boxs.f +
-</code> +
- +
- +
- +
-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 [[http://playground.metr.ou.edu/wiki/wiki.cgi/BrianFiedler|BrianFiedler]])+
  
  
  
  
tornadobox.1742332959.txt.gz · Last modified: by admin
Donate Powered by PHP Valid HTML5 Valid CSS Driven by DokuWiki