MATLAB VISUALIZATION TOOLS FOR CARPETcode

Gian Mario Manca
http://www.fis.unipr.it/~gianmario.manca
August 9th 2006

These scripts are distributed under the GNU General Public License (GPL). They might be released under the GNU Lesser General Public License (LGPL) in the future, to match the distribution terms of Cactus.

These scripts have to be considered just a starting point to facilitate the analysis and the representation of the 3D output of simulations made using Carpet.

Please use them, ask me about them, modify them, let the community (or at least me) know of your new scripts.

Just a little video (2.2Mb 1200x900) of a barmode instability.
It represents the mass density at z=0 and was made using CARPETsurfZ. With a little modification to add a light and emphasize the matter leaving the grid. Also the color-scale was fixed using the caxis matlab command. The simulation ran on 24 processors and the video was made starting from the 48 components of the 2 refinement level stored in the 24 chunked hdf5 files generated. PS: The title in the video is misleading: actually also the second refiniment level was used to generate the surface. By Gian Mario Manca.


Just a little video (3.7Mb 1200x900) of the first test of a 3D disk simulation.
It represents just a test showing the bounce of an unstable disk against the centrifugal barrier.
Improvement of the initial conditions and tests on the resolution effects are under way.
It was made using CARPETcontourfZ and you can see the mass density at z=0.
The color-scale was fixed using the caxis matlab command.
Simulation by Giovanni Corvino (mail: corvino at fis.unipr.it).


How did you make the movies?
Look also if you find something interesting here: Matlab file exchange.
Just to give you an example: try the Interactive Colorbar in combination with CARPETpcolorZ.m.

Principal scripts (download)


Tar.gz containing all the scripts:
CARPET_HDF5_MATLAB_scripts_12_03_2007.tar.gz

Older versions:
CARPET_HDF5_MATLAB_scripts_11_10_2006.tar.gz
CARPET_HDF5_MATLAB_scripts_9_8_2006.tar.gz


CARPETListDatasetHDF5.m
CARPETLoadDatasetHDF5.m
CARPETinterp3.m
CARPETcontourfZ.m
CARPETcontourfY.m
CARPETcontourfX.m
CARPETpcolorZ.m
CARPETpcolorX.m
CARPETpcolorY.m
CARPETisosurface.m

Variants (You can create a lot of variants, use different graphical commands, change the kind of the surfaces, insert lights, etc)

CARPETcontourfZ.m
CARPETmeshZ.m
CARPETsurfZ.m

Standard use

GO TO EXAMPLES FOR CHUNKED FILES
GO TO EXAMPLES FOR UNCHUNKED FILES

>> str=CARPETListDatasetHDF5('data','phi',1)

It creates a structure with all the information about the unchunked (chuncked) HDF5 file phi.h5 (phi.file_0.h5 phi.file_1.h5 ...) which is in the directory data. You have to use it ONLY ONCE before starting loading data.
  If you want you can also save the structure in a file typing ">> save filename.mat str" and reload it typing ">> load filename.mat" when you reopen your matlab for a new working session.
  Use an absolute path for the directory parameter if you want be able to load data from a directory different from that where the matlab was when you created the structure.

Remember that the script uses the name string (in the example "\WAVETOY::phi it=80 tl=0 rl=0 c=1") which characterizes every dataset. If you have problems check if in your hdf5 file (using for example the h5ls command of the hdf5 library) the it (iteration), rl (ref. level) and c (component) are in a different order or the string contains other information that could alter the result of the sscanf command at the interior of the script.

>> data=CARPETLoadDatasetHDF5(str,5)

It loads (USING the structure created by CARPETListDatasetHDF5) the 3D data of all the components of each refinement level of a specific iteration and put it in a structure. (READ AN IMPORTANT NOTE BELOW IF YOU HAVE PROBLEMS)

>> CARPETpcolorZ(data,valz)

It creates (USING the structure created by CARPETLoadDatasetHDF5) a graphical representation of the plan with z=valz coordinate.

If you have a component with a very high resolution it could appear black. Try to zoom on it or try shading interp or shading flat instead of shading faceted at the interior of the script. (The same for for CARPETsurf and CARPETmesh)

Try to change the extrema of the caxis command if you want to change the colorscale and to focus on a particular interval of values.

>> CARPETisosurface(data,var_value,transparency)

It reconstructs the surfaces where the quantity stored in data has the same value and superimposes the result obtained for each component in a matlab figure.

>> interpolated_values=CARPETinterp3(data,vector_x,vector_y,vector_z);

It creates a (nx*ny*nz)-matrix filled with interpolated values using the best ref. level for each point. (nx=length of vector_x ...)(You can change at the interior of the script the parameters of the interp3 function)

>> help CARPETListDatasetHDF5

REMEMBER: LOOK AT THE HELP!

   [str]= CARPETListDataSetsHDF5(direc,var,unc)
  
   EXAMPLE PARAMETERS
   direc= 'data/hdf5' (if you have your hdf5 file sin data/hdf5)
   var=   'rho'  -> look for the file rho.h5 in data/hdf5
                    for unchunked files
                 -> look for the file rho.file_0.h5 rho.file_1.h5 etc
                    in data/hdf5 for chuncked files
   unc=   0 -> chunked;  1-> unchunked
 
eems   CREATE a list of all the datasets in a given HDF5 file
 
    str.names   = ALL THE DATASETS' NAMES
    str.rl      = (vector) number of refinement levels for each iteration
    str.c       = (matrix) number of components for each iteration and ref-level
    str.it      = (vector) iterations EX: [0 20 40 60]
    str.index   = (vector) index corresponding to an iteration EX: [1 2 3 4]
    str.num2name= (matrix) gives the index of the name in str.names
                  corresponding to an iteration, ref. level, and component 
    str.ranges  = [first iterations, last iterations, step];
                 
 
    str.filename    = Filename;
    str.info        = info from hdf5info;
 
  HOW TO USE
  for a=str.index
      for b=1:str.rl(a)
          for c=1:str.c(a,b)
              disp(str.names{str.num2name(a,b,c)})
          end
      end
  end
 
   Author: Gian Mario Manca, manca@fis.unipr.it
           http://www.fis.unipr.it/~manca

EXAMPLES

UNCHUNKED FILES

>> str=CARPETListDatasetHDF5('data','phi',1)

Found data/phi.h5

str =

        info: {1x156 cell}
    num2file: [1x156 double]
    filename: {'data/phi.h5'}
          rl: [26x1 double]
           c: [26x3 double]
    num2name: [26x3x3 double]
          it: [1x26 double]
       index: [1x26 double]
        time: [1x26 double]
    timeMSEC: [1x26 double]
      ranges: [0 500 20]
       names: {1x156 cell}

>> data=CARPETLoadDatasetHDF5(str,5)

TIME (ms) = 0.036941 ms & TIME (C.U.) = 7.5 C.U.
/WAVETOY::phi it=80 tl=0 rl=0 ==> (61x61x61)  (1x1x1) X(-30,30) Y(-30,30) Z(-30,30)
/WAVETOY::phi it=80 tl=0 rl=1 c=0 ==> (21x31x21)  (0.5x0.5x0.5) X(-10,0) Y(-10,5) Z(-5,5)
/WAVETOY::phi it=80 tl=0 rl=1 c=1 ==> (20x21x21)  (0.5x0.5x0.5) X(0.5,10) Y(0,10) Z(-5,5)
/WAVETOY::phi it=80 tl=0 rl=1 c=2 ==> (6x7x9)  (0.5x0.5x0.5) X(0.5,3) Y(-13,-10) Z(-2,2)
/WAVETOY::phi it=80 tl=0 rl=2 c=0 ==> (7x9x17)  (0.25x0.25x0.25) X(0.5,2) Y(2,4) Z(-2,2)
/WAVETOY::phi it=80 tl=0 rl=2 c=1 ==> (10x8x17)  (0.25x0.25x0.25) X(0.75,3) Y(4.25,6) Z(-2,2)

data =

          rl: 3
           c: [1 3 2]
    num2name: [3x3 double]
        time: 7.5000
    timeMSEC: 0.0369
        xmin: -30
        ymin: -30
        zmin: -30
        xmax: 30
        ymax: 30
        zmax: 30
        vmax: 0.0685
        vmin: -0.1308
       names: {1x6 cell}
          dt: [1x6 struct]
        info: {1x6 cell}

IMPORTANT: If you have problems with the dimensions (X and Y dimension exchanged) look at the interior of the script and try to uncomment the line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% d.dt(num).val=permute(d.dt(num).val,[2 1 3]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
It depends on the hdf5 library and on the platform the way the data are loaded in a variable by the hdf5read matlab command.

Look also at the lines where the script extracts the information about deltaX, deltaY, deltaZ, and the coordinates of the left-bottom corner of the component:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check these lines if you have problems.
% Sometimes they are attributes 9 and 10, sometimes 1 and 2, sometimes who knows...
% It depends on your carpet version.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
left_bottom =str.info{str.num2name(index,i,j)}(9).Value; %d.attr(9).Value;
delta =str.info{str.num2name(index,i,j)}(10).Value; %d.attr(10).Value;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>> CARPETpcolorZ(data,0)

pcolor_zoom.jpg

>> CARPETcontourfZ(data,0)

contourf.jpg

>> CARPETisosurface(data,-0.03,0.2)

To change the colors of the surfaces in different components of a ref. level (as in this example) modify them using the menus of the matlab figures or change the script if you want something automatic.

isosurface.jpg

>> CARPETsurfZ(data,0)

surf.jpg surf2.jpg

CHUNKED FILES

>> str=CARPETListDatasetHDF5('data','phi',0)

Found data/phi.file_0.h5
Found data/phi.file_1.h5
Found data/phi.file_2.h5
Found data/phi.file_3.h5

str =

        info: {1x48 cell}
    num2file: [1x48 double]
    filename: {1x4 cell}
          rl: [4x1 double]
           c: [4x3 double]
    num2name: [4x3x4 double]
          it: [0 20 40 60]
       index: [1 2 3 4]
        time: [0 1.8750 3.7500 5.6250]
    timeMSEC: [0 0.0092 0.0185 0.0277]
      ranges: [0 60 20]
       names: {1x48 cell}

>> data=CARPETLoadDatasetHDF5(str,3)

TIME (ms) = 0.018471 ms & TIME (C.U.) = 3.75 C.U.
/WAVETOY::phi it=40 tl=0 rl=0 c=0 ==> (32x61x32)  (1x1x1) X(-30,1) Y(-30,30) Z(-30,1)
/WAVETOY::phi it=40 tl=0 rl=0 c=1 ==> (32x61x31)  (1x1x1) X(-30,1) Y(-30,30) Z(0,30)
/WAVETOY::phi it=40 tl=0 rl=0 c=2 ==> (31x61x32)  (1x1x1) X(0,30) Y(-30,30) Z(-30,1)
/WAVETOY::phi it=40 tl=0 rl=0 c=3 ==> (31x61x31)  (1x1x1) X(0,30) Y(-30,30) Z(0,30)
/WAVETOY::phi it=40 tl=0 rl=1 c=0 ==> (23x16x13)  (0.5x0.5x0.5) X(-10.5,0.5) Y(-10.5,-3) Z(-10.5,-4.5)
/WAVETOY::phi it=40 tl=0 rl=1 c=1 ==> (23x16x12)  (0.5x0.5x0.5) X(-10.5,0.5) Y(-10.5,-3) Z(-5,0.5)
/WAVETOY::phi it=40 tl=0 rl=1 c=2 ==> (23x9x23)  (0.5x0.5x0.5) X(-10.5,0.5) Y(-3.5,0.5) Z(-10.5,0.5)
/WAVETOY::phi it=40 tl=0 rl=1 c=3 ==> (13x13x13)  (0.5x0.5x0.5) X(4.5,10.5) Y(4.5,10.5) Z(4.5,10.5)
/WAVETOY::phi it=40 tl=0 rl=2 c=0 ==> (23x16x13)  (0.25x0.25x0.25) X(-5.25,0.25) Y(-5.25,-1.5) Z(-5.25,-2.25)
/WAVETOY::phi it=40 tl=0 rl=2 c=1 ==> (23x16x12)  (0.25x0.25x0.25) X(-5.25,0.25) Y(-5.25,-1.5) Z(-2.5,0.25)
/WAVETOY::phi it=40 tl=0 rl=2 c=2 ==> (23x9x23)  (0.25x0.25x0.25) X(-5.25,0.25) Y(-1.75,0.25) Z(-5.25,0.25)
/WAVETOY::phi it=40 tl=0 rl=2 c=3 ==> (11x11x11)  (0.25x0.25x0.25) X(5.75,8.25) Y(5.75,8.25) Z(5.75,8.25)

data =

          rl: 3
           c: [4 4 4]
    num2name: [3x4 double]
        time: 3.7500
    timeMSEC: 0.0185
        xmin: -30
        ymin: -30
        zmin: -30
        xmax: 30
        ymax: 30
        zmax: 30
        vmax: 0.1096
        vmin: -0.3991
       names: {1x12 cell}
          dt: [1x12 struct]
        info: {1x12 cell}

>> CARPETpcolorZ(data,0)

pcolor_ck.jpg

>> CARPETsurfZ(data,0)

surf2_ck.jpg

INTERPOLATION USING CARPETinterp3

>> x=[-10:0.2:10];
>> y=[-10:0.2:10];
>> z=[-1:1:1];
>> inter=CARPETinterp3(data,x,y,z);
>> surf(x,y,inter(:,:,1))
interp.jpg
>>[xg,yg,zg]=ndgrid(x,y,z);
>> temp=xg.*inter;
>> surf(x,y,temp(:,:,1))
interp2.jpg

INTEGRATION USING CARPETinterp3

(not at all efficient but it could be useful)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%str=CARPETListDatasetHDF5('data','phi',0)
%data=CARPETLoadDatasetHDF5(str,3)

dx=0.25;

x=[0:dx:10-dx];
y=[0:dx:10-dx];
z=[0:dx:10-dx];

total=0;

for i=[-30:10:20]
    for j=[-30:10:20]
        for k=[-30:10:20]
         %disp([num2str(i) ' ' num2str(j)  ' ' num2str(k)]);
         inter=CARPETinterp3(data,i+x,j+y,k+z);
         total=total+sum(sum(sum(inter)))*dx*dx*dx;
        end
    end
end

disp(total)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

With dx=0.5 total 115.947
With dx=0.25 total 115.974