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 ui,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: 2016/12/02 00:37:32 $, $Author: brandenb $, $Revision: 1.4 $