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/03 17:13] (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]]
Line 206: Line 208:
 $ ^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'  
 +'bf10',015000,'bF10'
 </code> </code>
 +$ 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 249:
   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: 2025/03/18 16:22 by admin
Donate Powered by PHP Valid HTML5 Valid CSS Driven by DokuWiki