Computer Physics: Exercise 13

Section 13 Exercise: Project 2, Part 2

Today's exercise is the second part of Project 2, the main theme of which is Visualization. In this part, you will learn hot to visualize and study snapshots of the same type of experiment that you started preparing in the previous exercise. You are going to do this using Paraview and IDL. Please note that there are no questions (and hence no credits) for this exercise, but since you need the techniques taught here to produce material for the report (which gives a lot of credit!), it is important to work also on this exercise.

If you haven't yet started the magnetic relaxsation experiment at DCSC-KU, start that tasks first, and then come back and work on the exercises below, using example data available on DCSC-KU at ~kg/scratch/cf10/12/ until you have your own data ready.

Subsections	approx time
CVS and data access	1 min
Data conversion	10 min
Paraview data visualization	2+ hours
IDL	30 min
Home Work	1 hour

Notes:

It is time for the course evaluations. It is done through the Absolon page for the course and it is open until Sunday, June 20th. !

The Thursday night dead-lines do not apply to Project 2, and there are no specific internal dead-lines for the various parts. The dead-line for finishing the whole project is Friday, June 25th.

CVS and data access

[about 1 minute]

To check out this weeks material, do the following:

    ssh -Y CF_user@fend03.dcsc.ku.dk
    cd ~/users/"computer-physics-id"/ComputerPhysics
    cvs update -AdP

This will create the directory 13_Visualization where you will do most of the work during this part of the project. For the following parts of the exercise to work you most work on DCSC-KU. Then do

    cd 13_Visualization
    \rm -f data                                  # if your dataset is ready
    ln -sf ~/scratch/"computer-physics-id" data                   
    ln -sf ../12_Relaxation/Idl/*.pro ./

    \rm -f data                                  # if your dataset is not ready yet
    ln -sf ~kg/scratch/cf_12/ data      
    ls -l data/snapshot.dat
    ls -l data/

to have access to the subdirectory containing the data you have to visualize. The ls commands are there as a check that you really have access to the data set via the name data/snapshot.dat.

If you start by using the example data, then when you have your own dataset ready to work on just redo the commands above, linking "data" to your own directory ~/scratch/"computer-physics-id" instead of to the example data directory.

VNC setup

When working from your own laptop, or if you want to do this from home, you need to access one of the graphics nodes on the cluster to be able to run paraview with full utilisation of the graphics card. First you have to start a VNC session on a graphical node. To do this first logon to fend03. Then log on to a graphical node:

  ssh $USER@fend03.dcsc.ku.dk
  llqstat | grep GPU           # we have reserved a few nodes for interactive usage, with job name GPU####
  ssh -p #### localhost        # where #### is the number indicated above

where the second line picks up the number of the nodes that are available for interactive usage. The information looks like this:

 mgnt03.2418724. troels   astro_gpu    GPU2222      1   8 240:0  R 06:16 node617

where 2222 in GPU2222 represents the port number after the -p option in the second ssh command above. The last entry node617 is the actual "name" of the machine you logon to. There will be more than one GPU#### node available. Each has 8 CPUs and 4 graphics cards. So please scatter randomly on these machines. Now start a VNC server on this node:

vncserver -geometry 1024x768 -depth 24 -localhost     # you may choose a smaller graphical window....

What you need to notice is that when starting the server provides a display number (in the form :N) that you need to remember for connecting (and later re-connecting) to the session with a vncviewer (see below). More information about setting up VNC is given under Course Material on Absalon.

Now exit from the GPU node and start a tunnel to the vncsession:

fend03>  ssh -p 2222 localhost -L 590#:localhost:590#         # replace 2222 with the port number and # with the display number

where # is the display number given by the initiation of the vncserver. To access the vnc session from your laptop, you need to setup a tunnel that connects to fend03 and then start the vncviewer:

laptop>  ssh -NfC CF_user@fend03.dcsc.ku.dk -L590#:localhost:590#
laptop>  vncviewer localhost:#

(or correspondingly with Putty or another SSH client).

Since you are all using the same login to DCSC-KU there is also a common VNC passwd: !CF_VNC_PASS. Give the passwd and your view of the vncsession running on the GPU machine will start. The graphical interface you will encounter is a simple one as neither kde nor gnome is installed on the graphical nodes.

Data conversion

[about 10 minutes]

Before you can use your own data from the relaxation simulations for visualization with PARAVIEW the data needs to be transformed into a format that Paraview reads.

NOTE: If initially you are using the example data you don't need to (and cannot) do the data conversion. Come back and do that later, with your own dataset.

Initially the data is written from the simulation code in a simple block format, where the data file contains 8 variables for each snapshot and on the order of 100 snapshots. Each of the variables are represented by N_x x N_y x N_z grid points.

Paraview is able to read a huge number of different formats. We have chosen to use their new standart format, XML, which is a platform independent open format. You can read more about it here

To transform data from a series of snapshots into XML format you need to execute the following command in IDL from the data directory:

    IDL> openr,10,'data/snapshot.dat'
    IDL> nn=192
    IDL> a=assoc(10,fltarr(nn,nn,nn))
    IDL> .c stagger
    IDL> for i=0,100 do stagger2vti,a,i,'data/snapshot',[1.,1.,1.] & print,i

Note: this takes about 1 minute per snapshot and each snapshot requires roughly 1.5 GB. You can follow the advance by the increase in the printed number. As soon as the first file is written you can continue the exercise below.

The process creates new files called snapshot.#.vti/, where # represent the snapshot number. Each file is selfcontained, including all the information need to reconstruct the data set when read into Paraview. Each file contains a number of processed variables (moved to center grid position), density, magnetic field vector, current vector and the velocity vector.

Paraview Visualization

[about 2+ hours]

Data visualization is a vital part in investigating numerical experiments. In this exercise you are going to learn how to make five different types of visualizations using the magnetic relaxation data. These represents five different ways of viewing the data:

Isosurface of magnitude of magnetic field
A translucent volume rendering of the current distribution for a not too early snapshot
Magnetic field lines showing the magnetic structure in the vicinity of null points
Velocity vectors showing the flow pattern in the dynamical evolution
A 2D image showing cross section of total density, in a plane through one of the relaxation regions

Combining the information from these independent visualizations puts you in a position to discuss the underlying relaxation processes (we also provide a text page with additional background for the discussion).

One can say that "visualization is a way to convert on the order of 8*N³ floating point numbers into physical intuition of what's going on in the numerical experiment".

Let's get started on the process.

Start Paraview in the VNC session by typing (in an Xterm window):

  paraview &

To access a dataset to investigate, click on File and choose Open. A file browser appears, choose the *vti file you have just created in your data directory. You need to press Apply to actually read the data (takes some time).

Below we discuss five ways to investigate the data set. Here we only described the process broadly, while more detailed instructions are available through links to a TWiki page within each block.

Isosurface of magnetic energy density: Due to the trigonometic defintion of the magnetic field the field strength varies a lot throughout the domain. The places we are specially interested in are the locations where magnetic nulls are present. To find these we need to choose a low value of magnetic magnitude. The isosurface represent a 3D surface where the magnetic energy density B²/2 (and hence the field strength B) is constant. The low values show where the magnetic field is week and therefore possible locations of magnetic null points. To see how to include an Isosurface into the Paraview image check the notes under the section Isosurfaces.
Volume rendering of the current density: Some time into the simulation the current density changes significantly as the initial Lorentz forces collapse the magnetic field structure in local regions into flat structures. Volume rendering is a technique used for investigating the spatial distribution of a scalar quantity inside a given part of the computational domain. The complicated magnetic field configurations gives rise to the formation of multiple locations with high current density. To understand the underlying physics behind the process creating these you first need to isolate one of these and investigate it in more detail. Volume rendering assigns a specific color value to a given magnitude of the scalar field, By changing the range used and the transparency of the colors one obtains a possibility of looking through the data set and to identify regions with different density values. To get an interesting sample of the field quantity some playing around with the dynamical range of the color table and the transparency level is required. Details on how to do this is found here under the section Volume visualization.
Velocity field vectors: One thing is to see the current density distribution or collapsed isosurfaces of magnetic field magnitude, another is to figure out why it looks as it does. One way to collapse a magnetic field into local sheet like structures, is when the Lorentz force pushes the field together, forming a local compression region. Another is when the field is exposed to a large shear across a given surface. This can be obtained by having twisting motions across the fan plane in null points, or by shearing field lines over a very short distance. Which of these that are responsible may be investigated looking at the velocity flow in the vicinity of the collapsed isosurface of the magnetic magnitude. Different regions have different reasons for collapsing, so looking at the symmetries of the current concentration and the local velocity flow can tell you why/how the structure is formed. To see how to include velocity vectors into the Paraview image check the notes under the section Vector field arrows visualization.
Magnetic field line traces: Knowing where the magnetic field is weak in space, we want to see how the collapse has modified the magnetic field structure. As you may remember, the induction equation tells how the magnetic field is mainly advected with the flow velocities. As the collpase requires some type of significant distorsion of the magnetic field by advecting the field, the structure of the field will have changed significantly in these regions. So, what does the magnetic field look like in the vicinity of a null point as it starts collapsing in time? By tracing field lines from a small neighborhood of such a region you will be able to find out. Again details on doing this are given here, under the section Vector field lines visualization.
Image slice: One still unused quantity stored in the data file represents the mass density. One can either look at this using volume rendering, an isosurface, or by showing the density in a particular plane passing through one of the null points. For more help look at the hints given here, under the section Image slices.

You now have the ingredients necessary to investigate and discuss how the physics generates the restructuring of the magnetic field, removing lot of the complexity with time. If you have, until now, been mostly busy doing the practical work that resulted in the Paraview renderings, then take a more detailed look at the visualization sequence again, and let the selection of variables and their special visualization methods tell the story about how the relaxation of the plasma takes place. You need this understanding for writing the final report, so talk to the teachers if you are in doubt. There is also a page with background information and discussion of the physical principles at work.

Before exiting Paraview you should save the Paraview settings (parameters, colors, perspectives, etc). To do this choose File and in the submenu choose Save state. This opens a new window, where you can choose a file name (ending with pvsm) and a directory to save the file in. This session can then be restored at a later time, to reconstruct the visualisation you have just created. We will return to this in the final exercise.

Visualization and data analysis using IDL

[about 30 minites]

As you have seen earlier, IDL can also handle various forms of 3D visualization. Here we are not going to pursue this, but instead show how one can use another simple tool to look at animations of the data and save the result as a MPEG movie. Let's see how the magnetic energy density map looks in a different projection. To do this we need to do the following:

   IDL>  .run stagger
   IDL>  n=256                                            ; grid resolution
   IDL>  n_snap=???                                       ; snapshot number to look at [0,100]
   IDL>  close,1                                          ; close file if open
   IDL>  openr,1,'data/snapshot.dat'                      ; open the data file for reading
   IDL>  a=assoc(1,fltarr(n,n,n))                         ; associates the variable "a" with the file
   IDL>  tmp=sqrt(xup(a[8*(n_snap-1)+5])^2+yup(a[8*(n_snap-1)+6])^2+zup(a[8*(n_snap-1)+7])^2) ;  magnetic magnitude snapshot n_snap
   IDL>  print,max(tmp,min=min),min                       ;  prints the max and min values

   IDL>  xinteranimate, set=[n,n,n], /showload            ;  allocate space for the animation
   IDL>  for i=0,n-1 do $                                 ;  loop over one index to show the content of the whole array
   IDL>     xinteranimate, frame=i, $                     ;  start loading data in to memory
   IDL>     image=bytscl(tmp[*,*,i],min=0,max=??)         ;  data for each animation frame 
   IDL>  xinteranimate, /keep_pixmaps                     ;  activates the animation tool

From the xinteranimate interface it is possible to view the animation, and choose the speed and direction of the animation. But, what is equally important is that this tool includes the ability to save the image sequence as an MPEG movie.

Checkout the documentation page on xinteranimate. It contains information about how to change the specific setting on for instance the quality of the MPEG movie using keywords to the first call to xinteranimate.

Similar animations can be made for the time evolution of other scalar quantities, showing how the local density, magnetic energy density, ... are evolving in time. Here's an example:

   IDL>  xinteranimate, set=[n,n,100?], /showload           ;  allocate space for the animation
   IDL>  for i=0,n_snap-1 do begin xinteranimate, $         ;  loop over one index to show the content of the array
   IDL>    frame=i,$                                        ;  frame number
   IDL>    image=bytscl(total(a[8*i+4],1),min=.1,max=5)     ;  density projection (along x) for each animation frame
   IDL>                                                     ;  bytescl data with in the range .1:5
   IDL>  xinteranimate, /keep_pixmaps                       ;  activates the animation tool

Other information may be found, by looking for example at distribution functions and peak values as functions of time. Here is an example with the peak value of the total density:

   IDL>  rho_max=fltarr(n_snap)
   IDL>  for i=0,n_snap-1 do begin rho_max[i]=max(a[8*i+4])
   IDL>  plot, rho_max, xtitle="??",ytitle="??"

To compute the root-mean-square of the current for a given snapshot do

   IDL>  Jx=ddyup(a[8*n_snap+7])-ddzup(a[8*n_snap+6])
   IDL>  Jy=...
   IDL>  Jz=...
   IDL>  Jrms=sqrt(total(Jx^2+Jy^2+Jz^2)/(3*n_elements(Jx)))

For magnetic relaxation experiments one may also look into the distribution function of the density, velocity amplitude, or magnetic energy in the experiment. For doing this one needs to use a function called histogram.

   IDL>  rho=a[8*n_snap]                                    ; get the density
   IDL>  h=histogram(alog10(rho),binsize=??,min=??,max=??)  ; compute histogram of log(rho)
   IDL>  plot,h,psym=10,/ylog                               ; plot the histogram

It is also possible to see how the different physical quantities change with time by making a 3D fourier analysis of ia iven quantity for different times. For this you need to use the special a routine called power3D.

   IDL>  Bx=a[8*n_snap+5]                                   ; get the x component of the B-filed
   IDL>  power3d,Bx                                         ; calculate and plot the 3D power spectrum

power3d has two keywords, wavenumbers and spectrum that can be used to get the spectra out, such that one can collect them into a 2D array that can then show the time evolution of the fourier transform.

To look for null points inside the domain you can use this setup

  IDL> bv=fltarr(n,n,n,3)                  ; define an array to contain the magnetic field vector
  IDL> bv(*,*,*,0)=xup(a[*n_snap+5])
  IDL> bv(*,*,*,1)=...
  IDL> bv(.....
  IDL> r_null=nullfinder(bv)               ; find the nulls in the domain -- returns the grid positions

Check the number of null point positions by determining the dimensions of the r_null array. How does this number change with time?