MATLAB VISUALIZATION TOOLS FOR CARPETcode
Gian Mario MancaThese 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.
Principal scripts (download)
|
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 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)
>> 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