Thoughts

This page shows a revisit to Rayleigh-Benard convection as a tool for studying dust devils, but with the computations done with numerical python. It is a substantial upgrade to my old dust devil page.

To reproduce these simulations, you may use my ipython notebook: 3D convection, periodic.

To reproduce the visualizations, you will need my companion python programs.

In the notebook, in cell number 28, set nexp to 91:

# choose an integer, see choices for nexp below
nexp = 91

And then, also in cell 28, you should set tend=40.

elif nexp == 91: # Fiedler and Kanak, 2001
    ray_over_rayc = 4096
    hftop = 0
    gamma = 0
    wcinc = .3
    outdir = 'R4096g0w4'
    xmax,ymax,zmax = 4,4,1
    #Nx,Ny,Nz= 129,129,33 # the original hi-res
    Nx,Ny,Nz = 65,65,17 # acceptable level of resolution
    kp = Nz//8
    kz = Nz//8
    neutral_init=True
    random_init=False
    anticipated_wmax=4.
    tend=30.
    dplot=.1

If you then select 'Run All Cells“ for the ipython notebook, the simulation will be complete in about 20 minutes (with my i7 chip). But the simulation for the images here was with the hi-resolution, or Nx,Ny,Nz= 129,129,33. With those specs, the simulation may need 5 hours (and 9 GB of disk space).

In cell 30, the buoyancy is initialized with

bi = .01*np.cos(4*np.pi*xb/xmax) + .01*np.cos(4*np.pi*yb/ymax) 
     + .00001*np.cos(6*np.pi*xb/xmax - 2*np.pi*yb/ymax)

I am not sure why that is apparently different from the idealized choices in Fiedler and Kanak, but it makes the point about symmetry breaking needed to produce the vertical vorticity of the extreme events, which are akin to dust devils.

Try

./hist.py R4096g0w4 1 400

if all went okay,

./hist.py R4096g0w4 1 40 > R4096g0w4.hist

0

And then

./plothist.py -s R4096g0w4.hist 

That should produce the following image, which shows a dust devil event just prior to t=30.

some extreme values at the surface

mkdir pngs
./npy2contour.py -o pngs R4096g0w4 000001 000400
./maketitle.py -o titleh.png
./png2mp4.py pngs/0 -t titleh.png

That will produce video.mp4, showing twin dust devils peaking at t=29.8:

The extreme event (a.k.a. dust devils) in 3-D, with isosurfaces.

mkdir vista
./npy2vista.py  R4096g0w4 000270 000320 -o vista
./maketitlevista.py
./png2mp4.py vista/0 -t titlevista.png

./npy2vista.py -s R4096g0w4 000298
open scene.html

Use your mouse, have fun …

$\zeta=-40$ (purple) , $p=-2$ (green), and $s$ at the surface

I did not finish implementing the command line options for npy2vista.py, but the source code can be easily edited to select the desired isosurfaces. (You may want to copy npy2vista.py to a file of some other name, and edit that…)

Change these lines:

#    plotter.add_mesh(contours_w_neg, color="lightblue", opacity=0.5, label="w = -0.5")
#    plotter.add_mesh(contours_w_pos, color="green", opacity=0.5, label="w = 0.5")
    plotter.add_mesh(contours_p1, color="green", opacity=0.5, label="p = -2")
#    plotter.add_mesh(contours_p2, color="purple", opacity=1.0, label="p = -3")
#    plotter.add_mesh(contours_p3, color="yellow", opacity=.5, label="s =  2")
#    plotter.add_mesh(contours_p4, color="red", opacity=.5, label="c =  6")
    plotter.add_mesh(contours_p5, color="purple", opacity=.5, label="$\zeta$ =  -40")
    plotter.add_mesh(edges, color='black', line_width=3, label="Domain Boundaries")
    plotter.add_mesh(origin_marker, color="yellow", label="Origin Marker")
    plotter.add_text(time,color='pink')
#    plotter.add_mesh(contour_2d, cmap="coolwarm", opacity=0.8, label="2D Contour (p at z=0)")
    plotter.add_mesh(slice_plane, scalars="s", cmap="OrRd", opacity=1.0, clim=[0,3.0 ])

to

    plotter.add_mesh(contours_w_neg, color="lightblue", opacity=0.5, label="w = -0.5")
    plotter.add_mesh(contours_w_pos, color="green", opacity=0.5, label="w = 0.5")
#    plotter.add_mesh(contours_p1, color="green", opacity=0.5, label="p = -2")
#    plotter.add_mesh(contours_p2, color="purple", opacity=1.0, label="p = -3")
#    plotter.add_mesh(contours_p3, color="yellow", opacity=.5, label="s =  2")
    plotter.add_mesh(contours_p4, color="red", opacity=.5, label="c =  6")
    plotter.add_mesh(contours_p5, color="purple", opacity=.5, label="$\zeta$ =  -40")
    plotter.add_mesh(edges, color='black', line_width=3, label="Domain Boundaries")
    plotter.add_mesh(origin_marker, color="yellow", label="Origin Marker")
    plotter.add_text(time,color='pink')
#    plotter.add_mesh(contour_2d, cmap="coolwarm", opacity=0.8, label="2D Contour (p at z=0)")
    plotter.add_mesh(slice_plane, scalars="s", cmap="OrRd", opacity=1.0, clim=[0,3.0 ])

And again:

./npy2vista.py -s R4096g0w4 000298
open scene.html

$w=.5$ (green), $w=-.5$ (grey) , $c=6$ (red)

The suitably renamed scene.html can then be uploaded to a website, as done here.