MATLAB VISUALIZATION TOOLS FOR CARPETcode

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.

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.mCARPETmeshZ.m CARPETsurfZ.m |

Standard use

GO TO EXAMPLES FOR CHUNKED FILESGO 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

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)

>> CARPETcontourfZ(data,0)

>> 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.

>> CARPETsurfZ(data,0)

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)

>> CARPETsurfZ(data,0)

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))

>>[xg,yg,zg]=ndgrid(x,y,z); >> temp=xg.*inter; >> surf(x,y,temp(:,:,1))

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