Instructions for data access

For people with cvs access; see details about the cvs repository. Do the following only once per computer:

export CVSROOT=:pserver:$USER@130.242.129.93:/var/cvs/brandenb
cvs login

To get the stuff, say, e.g.,

cvs co axel/forced/Bkm1Pot/KH1152tnuk4b_sig0_PrM100
cd !$
mv data data2
ln -s /scratch/axbr9098/pencil-code/axel/forced/Bkm1Pot/KH1152tnuk4b_sig0_PrM100 data
idl
pc_read_var,obj=var,/trimall,/bb,proc=11
help,var.bb

Point checks

To make sure you can read and process the data correctly, here some elementary tests. Start with the data from processor zero (data/proc0). The velocity array, for example, has the size

IDL> help,var.uu
    FLOAT     = Array[1152, 36, 32, 3]

The calculation of a few derived quantities is readily available via the /magic option, i.e.,

pc_read_var,/magic,/trimall,obj=var,variables=['uu','bb','uij','aa'],proc=0

But let's start verifying that the velocity is ok. The first few data points are:

IDL> print,var.uu[0:4,0:5,0,0] 
  3.51875e-05  3.70896e-05  6.98211e-05 -2.07822e-05  3.59904e-05
  0.000241414  0.000211438  0.000178087  0.000101530  5.43811e-06
  0.000326302  0.000298810  0.000272592  0.000217128  0.000127044
  0.000336363  0.000324796  0.000301157  0.000274452  0.000205764
  0.000298453  0.000307917  0.000308744  0.000286057  0.000237210
  0.000230666  0.000257294  0.000262968  0.000264212  0.000217930

And now the same thing for the magnetic vector potential:

IDL> print,var.aa[0:4,0:5,0,0]
  0.000647369  0.000537547  0.000431973  0.000354911  0.000336831
  0.000787110  0.000693840  0.000591880  0.000498778  0.000454261
  0.000898826  0.000820049  0.000726497  0.000633220  0.000572499
  0.000961314  0.000897281  0.000819807  0.000743500  0.000692698
  0.000962660  0.000916380  0.000863078  0.000819884  0.000803676
  0.000905840  0.000880839  0.000858411  0.000859878  0.000895758

Next, the same for the magnetic field:

IDL> print,var.bb[0:4,0:5,0,0]
  -0.00212512  -0.00188570  -0.00170409  -0.00328576  -0.00354289
  -0.00157956  -0.00118189  -0.00108118  -0.00118743  -0.00260845
  -0.00125290  -0.00118169  -0.00106216  -0.00113669  -0.00128127
 -0.000737481 -0.000872292 -0.000830125 -0.000869349 -0.000968076
 -0.000446811 -0.000655431 -0.000733316 -0.000620224 -0.000677850
 -0.000305913 -0.000450134 -0.000528719 -0.000587892 -0.000583366

To get hold of the rate of strain matrix, you'd first need the velocity gradient matrix. With idl, it is quite simple. We have already computed uij, which is the same as u_i,j. It gives

IDL> print,var.uij[0:4,0:5,0,0,0]
   0.00213247   0.00567780  -0.00812415  -0.00286018  0.000661866
  -0.00126566  -0.00588929  -0.00905461   -0.0209631   0.00160756
  -0.00301168  -0.00496871  -0.00661058   -0.0136346   -0.0197034
 -0.000390647  -0.00378170  -0.00395457  -0.00812884   -0.0160216
   0.00392791  0.000736585  -0.00164115  -0.00618087   -0.0126964
   0.00610974   0.00259065   0.00116283  -0.00346733   -0.0122392

IDL> print,var.uij[0:4,0:5,0,1,2]
  -0.00575988   0.00497872   0.00935728  -0.00377012    0.0325513
  -0.00726452   0.00244051   0.00724671   0.00333783    0.0105773
  -0.00825590  -0.00498893   0.00421164  -0.00125827    0.0128633
  -0.00958705  -0.00612219  0.000870539   0.00107675   0.00707754
   -0.0105986  -0.00789873  -0.00249023   0.00141597   0.00816002
   -0.0110035  -0.00737753  -0.00562428  -0.00102747   0.00872247

Next, let me first collapse all points in x,y,z into a single string, i.e.,

s=size(var.uij) & nx=s[1] & ny=s[2] & nz=s[3]
uij=reform(var.uij,1L*nx*ny*nz,3,3) & sij=uij
for i=0,2 do begin & for j=0,2 do begin & sij(*,i,j)=.5*(uij(*,i,j)+uij(*,j,i))

Let's now check that we are still on the same page:

IDL> help,sij
SIJ             FLOAT     = Array[1327104, 3, 3]
IDL> print,reform(sij[0,*,*])
   0.00213247    0.0144861   0.00569494
    0.0144861    0.0156897   -0.0147725
   0.00569494   -0.0147725   -0.0170805

This is thus sij on the first meshpoint. It is symmetric, as it should be. And here another point:

IDL> print,reform(sij[17,*,*])
   -0.0134711 -0.000133410    0.0176575
 -0.000133410    0.0242518   -0.0324144
    0.0176575   -0.0324144   -0.0103086

Now the eigenvalues, which are hidden in the last row. Check first

vecs=eigvec3_arr(sij)

So the 3 eigenvalues are available via

vecs=eigvec3_arr(sij)
IDL> help,vecs
VECS            COMPLEX   = Array[1327104, 3, 4]
pdf,n=500,reform(vecs(*,0,3)),xe1,ye1
pdf,n=500,reform(vecs(*,1,3)),xe2,ye2
pdf,n=500,reform(vecs(*,2,3)),xe3,ye3

where I called my private pdf routine. If you plot the 3 histograms, you should get something like in one of our handouts. Here some explicit values for the three eigenvalues at the first and 17th points:

IDL> print,reform(float(vecs[0,*,3]))
    0.0271460  0.000545120   -0.0269494
IDL> print,reform(float(vecs[17,*,3]))
    0.0452250  -0.00573060   -0.0390223

So the largest eigenvalue comes first. The corresponding eigenvectors are

IDL> print,reform(float(vecs[0,*,0:2]))
     0.447715     0.863852    -0.230893
     0.814464    -0.287398     0.504035
     0.369053    -0.413717    -0.832249

IDL> print,reform(float(vecs[17,*,0:2]))
     0.162978    -0.828660     0.535500
     0.836826     0.403615     0.369889
     0.522648    -0.387837    -0.759224

Here, print,reform(float(vecs[0,*,0])) is the eigenvector belonging to the largest eigenvalue,

IDL> print,reform(float(vecs[0,*,0]))
     0.447715     0.863852    -0.230893

The one belonging to the intermediate eigenvalue is

IDL> print,reform(float(vecs[0,*,1]))
     0.814464    -0.287398     0.504035

etc. The 3 eigenvectors are (or course) orthogonal to each other. em compute the alignments with the magnetic B-vector, you could do:

bb=reform(var.bb,1L*nx*ny*nz,3)
b=1./sqrt(dot2(var.bb))
cosb1=b*dot(bb,vecs[*,*,0]) & pdf,fmin=0,fmax=1,n=500,abs(cosb1),xb1,yb1
cosb2=b*dot(bb,vecs[*,*,1]) & pdf,fmin=0,fmax=1,n=500,abs(cosb2),xb2,yb2
cosb3=b*dot(bb,vecs[*,*,2]) & pdf,fmin=0,fmax=1,n=500,abs(cosb3),xb3,yb3

If you plot this

plot,xb1,yb1,xr=[0,1],yr=[0,2.5]
oplot,xb2,yb2
oplot,xb3,yb3

You should get a familiar plot.

$Date: 2017/08/15 12:32:19 $, $Author: brandenb $, $Revision: 1.5 $