Plotting 2D column-shaped results with python

Ok, so you have your computational output of a set of 2D points. You have been lazy and done the obvious stuff: arrange them in columns, with the first one being the x coordinates, second one the y coordinate, then come the fields, which may be scalar (one column each), or vector (two columns each, one for each coordinate). How to visualize them?

ipython --pylab

Then, read the data

In : dt = loadtxt( ‘1/mesh.dat’ )

In : shape( dt )
Out: (1024, 27)

Notice the last command tells us we have 1024 data points, and 27 fields (well, 25 + positions). For convinience, assign columns to arrays:

In : x=dt[:,0]; y=dt[:,1]; al=dt[:,4]

Now x and y are positions, and “al” is the scalar field for the fifth column (number 4, since counters start at 0 in python).

To visualize the positions,

In : scatter( x , y )

A scalar field may be visualized with a color map:

In : scatter( x , y , c = al )

The “c=” means the color is taken from field al. One may fiddle with colormaps and symbol sizes:

In : scatter( x , y , c = al , cmap= plt.cm.Blues, s=8 )

To know the range we are plotting, produce a color bar:

In : colorbar()

Remember each plotting is overlaid on the previous one, so it is necessary to blank the plot from time to time:

In : clf()

For vector fields, assign coordinates to two separate arrays:

In : vx=dt[:,8]; vy=dt[:,9]

Then, use “quiver” to get a vector plot:

In : quiver( x, y, vx , vy )

A markdown test

This is a test of markdown blog writing. The writing comes straight from my website on CFD methods. These were written in markdown under reveal.js, for quick and nice lecture slides. Some changes had to be made:

• LaTeX must start as “dollar sign latex” … “dollar”
• Links to local files (such as pictures) don’t work
• Lists (such as this one) do not seem to work well

Muy a menudo, se parte de las EDPs, conocidas, por ejempo: $\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0$

Estas se discretizan: sustituyendo las derivadas por diferencias.

Sin embargo, este es un proceso de ida y vuelta, porque
las EDPs se deducen a nivel discreto.

Deducción

Se suponen cambios de un campo $u$ sólo
en la dirección $x$    El cambio en la cantidad total $A \Delta x \, u_i$ será: $\frac{d }{d t} (A \Delta x \, u_i ) = \Phi_{i-1/2} - \Phi_{i+1/2}$

Flujos, convección

Antes los flujos por las caras venían dados por: $\Phi_{i-1/2} = A c \, u_{i-1/2}$ $\Phi_{i+1/2} = A c \, u_{i+1/2}$ $\frac{d }{d t} (A \Delta x \, u_i ) = A c \, u_{i-1/2} - A c \, u_{i+1/2}$

Sparse Poisson problem in eigen Back to scientific computing. Lately, I have been using the Eigen libraries for linear algebra. They were recommended by the CGAL project, and they indeed share a common philosophy. Thanks to the rather high C++ level they can accomplish this sort of magic:

int n = 100;
VectorXd x(n), b(n);
SpMat A(n,n);

fill_A(A);

fill_b(b);

// solve Ax = b
Eigen::BiCGSTAB<SpMat> solver;

solver.compute(A);

x = solver.solveWithGuess(b,x0);

Notice that A is a sparse matrix! I am next describing how to use this in order to solve the 1D Poisson equation.

OpenFOAM cheatlist

A quick cheatsheet for OpenFOAM. In italics, things that are useful but not part of OpenFOAM proper. Interesting read: The OpenFOAM Technology Primer

Sites

Shortcuts to directories

(type “alias” to reveal these)

• run (go to own’s running directory)
• foam
• foamfv
• foam3rdParty (hit <tab> for these longish commands!)
• tut
• app
• sol
• util
• lib
• src

Environment variables

• echo $FOAM_ <tab> (directories) • echo$WM_ <tab>  (building, aka compiling, settings)

Mesh

• Generation
• Import / export
• foamMeshToFluent
• fluentMeshToFoam, etc …
• Operations
• refineHexMesh
• transformPoints
• makeAxialMesh
• collapseEdges
• autoPatch
• mirorMesh
• Properties
• checkMesh

Fields

• setFields
• topoSet
• patchAverage
• patchIntegrate
• vorticity
• yPlusRAS
• yPlusLES
• boxTurb
• applyBoundaryLayer
• R
• wallShearStress

Help

• foamHelp (e.g. foamHelp boudary -field U)

Postprocessing

• sample
• paraview
• probeLocations

Solvers

• icoFoam
• interFoam
• many, many others

Running

• foamJob
• decomposePar
• reconstructPar
• mpirun
• nohup

Building

• foamNew source App …
• doxygen doxyfile
• gdb
• valgrind
• wmake
• wclean
• aliases for settings: wm32, wm64, wmSP, wmDP, wmSET, wmUNSET

Sorting eivenvectors in octave

A brief note. Doing ee= eig(m); produces a vector with all eigenvalues of matrix m. These are unsorted, so ee=sort(eig(m)); produces a vector with a sorted list of all eigenvalues of matrix m.

If we want the eigenvectors, [vv,eed]=eig(m); produces a diagonal matrix eed whose elements are the same as ee, and are unsorted. Who to produce a sorted vector with the eigenvalues, and re-order the eigenvectors accordingly? Thus:

[vv,eed]=eig(m);

[ee,perm] = sort(diag(eed));

vv=vv(:,perm);

from gnuplot to python

Part of my moving to python for science.

I have lately been using gnuplot to monitor the progress and end results of my simulations.  These are 2D hydrodynamic simulations, which involve 2D scalar and vector fields.

2D plots

With gnuplot, I used to run stuff like

plot ‘8000/points.dat’ u 1:2

for a 2D plot of a scalar field. On column 1 I have the x values of the coordinates and on column 2 the y values, so this just plots the coordinates. The “8000” is the directory for time step 8000, where the data file is saved.

OK, now running ipython –pylab I have to run

ss=scatter(dt[:,0],dt[:,1],c=dt[:,11],s=50) 2D scatter plot with colors

Important things to notice:

• python starts numbering c-style: at 0! so, despite pylab having been designed in order to mimic matlab syntax, there is a clear departure here.
• it’s easy to assign a color to the points, in this case column 12 (11 for python) is a density field. The default color scale looks ok (it goes from blue to red).
• the size is set to 50 for a nice full-screen view.
• New plots will appear on top of this one. If this is not desired, close the window, or run pylab.clf()

Vector plots

This is quite easy in gnuplot:

plot ‘8000/points.dat’ u 1:2:($5/10):($6/10) w vec

On columns 5 and 6 I store the x and y components of the velocity. In this example they are reduced ten-fold in order to hace a decent plot (and this has to set by hand…).

Now we may run:

st=quiver(dt[:,0],dt[:,1],dt[:,4],dt[:,5],dt)

which does the same. A bonus is that the vector length is automatically calculated from the average vector length. This may the right solution, if not the option “scale=xx” may be used. The higher xx is, the shorter the arrows (??). The fifth entry is a color code, which does not look to good at present because, I think, they are taken as RGB values, unlike with scatter.

3D plots

This one’s harder. To get a scatter 3D plot of three columns I used to run

splot ‘8000/points.dat’ u 1:2:12

Now, it’s more like

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()

ax = fig.add_subplot(111, projection=’3d’)

ax.scatter(dt[:,0],dt[:,1],dt[:,11])

fig.show()

That produces a quite nice scatter plot with tick marks on the three axes.

xubuntu installation

Notes on the installation of xubuntu 13.10, Saucy Salamander, on my work desktop.

• Long story short: problems with openSUSE (suddenly not booting…), sick and tired of the new KDE and stuff I never use anyway but slow down my system (nepomuk, akonandi, virtuoso, plasma thingies). Tried debian, and was pretty happy; unfortunately, openFOAM was a bit of a pain to install, due to issues with the libc libraries (that I sort of worked out by compiling from scratch). On top of that, the graphic driver did not suppor dual monitors, with which I can live, but paraview would not work (missing GLX, that stuff). Propietary ATI drivers also did not work.
• On the upside, I fell in love with the light-weight Xfce windows manager, so I tried xubuntu next. It is just a port to Xfce of a standard ubuntu distribution.
• Dowloaded iso image from site. I just took the suggested version, perhaps should have tried the LTS version instead 😛
• Burned onto DVD drive. The iso image is just too big for a CD.
• Booted with no problem onto live system, installed carefully, respecting existing /home. Network was impossible to configure manually at this stage.
• Installed onto hard-drive, booted from hard drive, everything was fine, with all my files and previous settings running. Now the network can be manually set.
• Tons of updates are now suggested. Some are related to graphic drivers (let’s keep our fingers crossed). Yep, X11 still works after restarting.
• Ubuntu software center (USC). First thing: paraview. It works!
• Next: openFOAM. Not included in USC, only the freefoam fork (I’ll try it some day). From the official download page, everything works as explained for ubuntu 13.10, codename saucy. No need to compile, which is such a pain (at least for the poor computer… it takes hours!). It is convenient to start this way, since tons of other stuff I’ll need later (g++, boost, the works) is also installed. In order for paraFoam to work, cmake should be installed.
• CGAL libs. Some helpful comments here. I didn’t need any special libs after having installed all the other stuff. Only cmake-gui, because for some reason the working directory was set to “4.2” in the latest 4.3 release (I don’t know if this is my fault or theirs). Ah, had to install liblapack-dev.
• octave and gnuplot. These come in USC
• ipython, scipy. Also in USC
• LaTeX. Just install beamer from USC, and a whole set of LaTeX things gets installed (latex, pdflatex, bibtex,…). texlive-latex-extra for the exam class. texlive-publishers for a bunch of interesting classes. texlive-lang-spanish also.
• Cloud stuff: spideroak, dropbox (in USC)
• Internet: chromium (in USC)
• vim and emacs (in USC)
• The classic xfig. Don’t forget gsfonts-x11 in order fonts work, and restarting the X system.