Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F72690994
io.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Jul 16, 14:29
Size
38 KB
Mime Type
text/x-c
Expires
Thu, Jul 18, 14:29 (2 d)
Engine
blob
Format
Raw Data
Handle
19072608
Attached To
R195 SCITAS application benchmarks
io.c
View Options
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <errno.h>
#ifdef HAVE_HDF5
#include <hdf5.h>
#endif
#include "allvars.h"
#include "proto.h"
/*! \file io.c
* \brief Routines for producing a snapshot file on disk.
*/
static
int
n_type
[
6
];
static
long
long
ntot_type_all
[
6
];
/*! This function writes a snapshot of the particle distribution to one or
* several files using the selected file format. If NumFilesPerSnapshot>1,
* the snapshot is distributed onto several files, several of them can be
* written simultaneously (up to NumFilesWrittenInParallel). Each file
* contains data from a group of processors.
*/
void
savepositions
(
int
num
)
{
double
t0
,
t1
;
char
buf
[
500
];
int
i
,
j
,
*
temp
,
n
,
filenr
,
gr
,
ngroups
,
masterTask
,
lastTask
;
t0
=
second
();
if
(
ThisTask
==
0
)
printf
(
"
\n
writing snapshot file...
\n
"
);
#if defined(SFR) || defined(BLACK_HOLES)
rearrange_particle_sequence
();
/* ensures that new tree will be constructed */
All
.
NumForcesSinceLastDomainDecomp
=
1
+
All
.
TreeDomainUpdateFrequency
*
All
.
TotNumPart
;
#endif
if
(
NTask
<
All
.
NumFilesPerSnapshot
)
{
if
(
ThisTask
==
0
)
printf
(
"Fatal error.
\n
Number of processors must be larger or equal than All.NumFilesPerSnapshot.
\n
"
);
endrun
(
0
);
}
if
(
All
.
SnapFormat
<
1
||
All
.
SnapFormat
>
3
)
{
if
(
ThisTask
==
0
)
printf
(
"Unsupported File-Format
\n
"
);
endrun
(
0
);
}
#ifndef HAVE_HDF5
if
(
All
.
SnapFormat
==
3
)
{
if
(
ThisTask
==
0
)
printf
(
"Code wasn't compiled with HDF5 support enabled!
\n
"
);
endrun
(
0
);
}
#endif
/* determine global and local particle numbers */
for
(
n
=
0
;
n
<
6
;
n
++
)
n_type
[
n
]
=
0
;
for
(
n
=
0
;
n
<
NumPart
;
n
++
)
n_type
[
P
[
n
].
Type
]
++
;
#ifdef CHECK_TYPE_DURING_IO
if
(
n_type
[
0
]
!=
N_gas
)
{
printf
(
"(%d) n_type[0]=%d N_gas=%d !!!
\n\n
"
,
ThisTask
,
n_type
[
0
],
N_gas
);
endrun
(
111000
);
}
#ifdef SFR
if
(
n_type
[
ST
]
!=
N_stars
)
{
printf
(
"(%d) n_type[ST]=%d N_gas=%d !!!
\n\n
"
,
ThisTask
,
n_type
[
ST
],
N_stars
);
endrun
(
111001
);
}
#endif
#endif
/* because ntot_type_all[] is of type `long long', we cannot do a simple
* MPI_Allreduce() to sum the total particle numbers
*/
temp
=
malloc
(
NTask
*
6
*
sizeof
(
int
));
MPI_Allgather
(
n_type
,
6
,
MPI_INT
,
temp
,
6
,
MPI_INT
,
MPI_COMM_WORLD
);
for
(
i
=
0
;
i
<
6
;
i
++
)
{
ntot_type_all
[
i
]
=
0
;
for
(
j
=
0
;
j
<
NTask
;
j
++
)
ntot_type_all
[
i
]
+=
temp
[
j
*
6
+
i
];
}
free
(
temp
);
/* assign processors to output files */
distribute_file
(
All
.
NumFilesPerSnapshot
,
0
,
0
,
NTask
-
1
,
&
filenr
,
&
masterTask
,
&
lastTask
);
fill_Tab_IO_Labels
();
if
(
All
.
NumFilesPerSnapshot
>
1
)
sprintf
(
buf
,
"%s%s_%04d.%d"
,
All
.
OutputDir
,
All
.
SnapshotFileBase
,
num
,
filenr
);
else
sprintf
(
buf
,
"%s%s_%04d"
,
All
.
OutputDir
,
All
.
SnapshotFileBase
,
num
);
ngroups
=
All
.
NumFilesPerSnapshot
/
All
.
NumFilesWrittenInParallel
;
if
((
All
.
NumFilesPerSnapshot
%
All
.
NumFilesWrittenInParallel
))
ngroups
++
;
for
(
gr
=
0
;
gr
<
ngroups
;
gr
++
)
{
if
((
filenr
/
All
.
NumFilesWrittenInParallel
)
==
gr
)
/* ok, it's this processor's turn */
write_file
(
buf
,
masterTask
,
lastTask
);
MPI_Barrier
(
MPI_COMM_WORLD
);
}
if
(
ThisTask
==
0
)
printf
(
"done with snapshot.
\n
"
);
t1
=
second
();
All
.
CPU_Snapshot
+=
timediff
(
t0
,
t1
);
}
/*! This function fills the write buffer with particle data. New output blocks
* can in principle be added here.
*/
void
fill_write_buffer
(
enum
iofields
blocknr
,
int
*
startindex
,
int
pc
,
int
type
)
{
int
n
,
k
,
pindex
;
float
*
fp
;
#ifdef LONGIDS
long
long
*
ip
;
#else
int
*
ip
;
#endif
#ifdef PERIODIC
FLOAT
boxSize
;
#endif
#ifdef PMGRID
double
dt_gravkick_pm
=
0
;
#endif
double
dt_gravkick
,
dt_hydrokick
,
a3inv
=
1
,
fac1
,
fac2
;
if
(
All
.
ComovingIntegrationOn
)
{
a3inv
=
1
/
(
All
.
Time
*
All
.
Time
*
All
.
Time
);
fac1
=
1
/
(
All
.
Time
*
All
.
Time
);
fac2
=
1
/
pow
(
All
.
Time
,
3
*
GAMMA
-
2
);
}
else
a3inv
=
fac1
=
fac2
=
1
;
#ifdef PMGRID
if
(
All
.
ComovingIntegrationOn
)
dt_gravkick_pm
=
get_gravkick_factor
(
All
.
PM_Ti_begstep
,
All
.
Ti_Current
)
-
get_gravkick_factor
(
All
.
PM_Ti_begstep
,
(
All
.
PM_Ti_begstep
+
All
.
PM_Ti_endstep
)
/
2
);
else
dt_gravkick_pm
=
(
All
.
Ti_Current
-
(
All
.
PM_Ti_begstep
+
All
.
PM_Ti_endstep
)
/
2
)
*
All
.
Timebase_interval
;
#endif
fp
=
CommBuffer
;
ip
=
CommBuffer
;
pindex
=
*
startindex
;
switch
(
blocknr
)
{
case
IO_POS
:
/* positions */
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
for
(
k
=
0
;
k
<
3
;
k
++
)
{
fp
[
k
]
=
P
[
pindex
].
Pos
[
k
];
#ifdef PERIODIC
boxSize
=
All
.
BoxSize
;
#ifdef LONG_X
if
(
k
==
0
)
boxSize
=
All
.
BoxSize
*
LONG_X
;
#endif
#ifdef LONG_Y
if
(
k
==
1
)
boxSize
=
All
.
BoxSize
*
LONG_Y
;
#endif
#ifdef LONG_Z
if
(
k
==
2
)
boxSize
=
All
.
BoxSize
*
LONG_Z
;
#endif
while
(
fp
[
k
]
<
0
)
fp
[
k
]
+=
boxSize
;
while
(
fp
[
k
]
>=
boxSize
)
fp
[
k
]
-=
boxSize
;
#endif
}
n
++
;
fp
+=
3
;
}
break
;
case
IO_VEL
:
/* velocities */
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
if
(
All
.
ComovingIntegrationOn
)
{
dt_gravkick
=
get_gravkick_factor
(
P
[
pindex
].
Ti_begstep
,
All
.
Ti_Current
)
-
get_gravkick_factor
(
P
[
pindex
].
Ti_begstep
,
(
P
[
pindex
].
Ti_begstep
+
P
[
pindex
].
Ti_endstep
)
/
2
);
dt_hydrokick
=
get_hydrokick_factor
(
P
[
pindex
].
Ti_begstep
,
All
.
Ti_Current
)
-
get_hydrokick_factor
(
P
[
pindex
].
Ti_begstep
,
(
P
[
pindex
].
Ti_begstep
+
P
[
pindex
].
Ti_endstep
)
/
2
);
}
else
dt_gravkick
=
dt_hydrokick
=
(
All
.
Ti_Current
-
(
P
[
pindex
].
Ti_begstep
+
P
[
pindex
].
Ti_endstep
)
/
2
)
*
All
.
Timebase_interval
;
for
(
k
=
0
;
k
<
3
;
k
++
)
{
fp
[
k
]
=
P
[
pindex
].
Vel
[
k
]
+
P
[
pindex
].
GravAccel
[
k
]
*
dt_gravkick
;
if
(
P
[
pindex
].
Type
==
0
)
{
fp
[
k
]
+=
SphP
[
pindex
].
HydroAccel
[
k
]
*
dt_hydrokick
;
}
}
#ifdef PMGRID
for
(
k
=
0
;
k
<
3
;
k
++
)
fp
[
k
]
+=
P
[
pindex
].
GravPM
[
k
]
*
dt_gravkick_pm
;
#endif
for
(
k
=
0
;
k
<
3
;
k
++
)
fp
[
k
]
*=
sqrt
(
a3inv
);
n
++
;
fp
+=
3
;
}
break
;
case
IO_ID
:
/* particle ID */
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
ip
++
=
P
[
pindex
].
ID
;
n
++
;
}
break
;
case
IO_MASS
:
/* particle mass */
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
P
[
pindex
].
Mass
;
n
++
;
}
break
;
case
IO_U
:
/* internal energy */
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
#ifdef ISOTHERM_EQS
*
fp
++
=
SphP
[
pindex
].
Entropy
;
#else
#ifdef MULTIPHASE
switch
(
SphP
[
pindex
].
Phase
)
{
case
GAS_SPH
:
*
fp
++
=
dmax
(
All
.
MinEgySpec
,
SphP
[
pindex
].
Entropy
/
GAMMA_MINUS1
*
pow
(
SphP
[
pindex
].
Density
*
a3inv
,
GAMMA_MINUS1
));
break
;
case
GAS_STICKY
:
*
fp
++
=
SphP
[
pindex
].
Entropy
;
break
;
case
GAS_DARK
:
*
fp
++
=
-
SphP
[
pindex
].
Entropy
;
break
;
}
#else
*
fp
++
=
dmax
(
All
.
MinEgySpec
,
SphP
[
pindex
].
Entropy
/
GAMMA_MINUS1
*
pow
(
SphP
[
pindex
].
Density
*
a3inv
,
GAMMA_MINUS1
));
#endif
#endif
n
++
;
}
break
;
case
IO_RHO
:
/* density */
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
SphP
[
pindex
].
Density
;
n
++
;
}
break
;
case
IO_HSML
:
/* SPH smoothing length */
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
SphP
[
pindex
].
Hsml
;
n
++
;
}
break
;
/*********************************************/
/* here it is the end of the minimal output */
/*********************************************/
case
IO_STAR_FORMATIONTIME
:
/* stellar formation time */
#ifdef OUTPUTSTELLAR_PROP
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
StP
[
P
[
pindex
].
StPIdx
].
FormationTime
;
n
++
;
}
#endif
break
;
case
IO_INITIAL_MASS
:
/* stellar formation time */
#ifdef OUTPUTSTELLAR_PROP
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
StP
[
P
[
pindex
].
StPIdx
].
InitialMass
;
n
++
;
}
#endif
break
;
case
IO_STAR_IDPROJ
:
/* stellar projenitor */
#ifdef OUTPUTSTELLAR_PROP
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
ip
++
=
StP
[
P
[
pindex
].
StPIdx
].
IDProj
;
n
++
;
}
#endif
break
;
case
IO_STAR_RHO
:
/* gas density around a star */
#ifdef OUTPUTSTELLAR_PROP
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
StP
[
P
[
pindex
].
StPIdx
].
Density
;
n
++
;
}
#endif
break
;
case
IO_STAR_HSML
:
/* SPH smoothing length of a star */
#ifdef OUTPUTSTELLAR_PROP
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
StP
[
P
[
pindex
].
StPIdx
].
Hsml
;
n
++
;
}
#endif
break
;
case
IO_STAR_METALS
:
/* stars metal */
#ifdef OUTPUTSTELLAR_PROP
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
for
(
k
=
0
;
k
<
NELEMENTS
;
k
++
)
{
fp
[
k
]
=
StP
[
P
[
pindex
].
StPIdx
].
Metal
[
k
];
}
n
++
;
fp
+=
NELEMENTS
;
}
#endif
break
;
case
IO_METALS
:
/* gas metal */
#ifdef OUTPUTSTELLAR_PROP
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
for
(
k
=
0
;
k
<
NELEMENTS
;
k
++
)
{
fp
[
k
]
=
SphP
[
pindex
].
Metal
[
k
];
}
n
++
;
fp
+=
NELEMENTS
;
}
#endif
break
;
case
IO_POT
:
/* gravitational potential */
#ifdef OUTPUTPOTENTIAL
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
P
[
pindex
].
Potential
;
n
++
;
}
#endif
break
;
case
IO_ACCEL
:
/* acceleration */
#ifdef OUTPUTACCELERATION
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
for
(
k
=
0
;
k
<
3
;
k
++
)
fp
[
k
]
=
fac1
*
P
[
pindex
].
GravAccel
[
k
];
#ifdef PMGRID
for
(
k
=
0
;
k
<
3
;
k
++
)
fp
[
k
]
+=
fac1
*
P
[
pindex
].
GravPM
[
k
];
#endif
if
(
P
[
pindex
].
Type
==
0
)
for
(
k
=
0
;
k
<
3
;
k
++
)
{
fp
[
k
]
+=
fac2
*
SphP
[
pindex
].
HydroAccel
[
k
];
}
fp
+=
3
;
n
++
;
}
#endif
break
;
case
IO_DTENTR
:
/* rate of change of entropy */
#ifdef OUTPUTCHANGEOFENTROPY
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
SphP
[
pindex
].
DtEntropy
;
n
++
;
}
#endif
break
;
case
IO_TSTP
:
/* timestep */
#ifdef OUTPUTTIMESTEP
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
(
P
[
pindex
].
Ti_endstep
-
P
[
pindex
].
Ti_begstep
)
*
All
.
Timebase_interval
;
n
++
;
}
#endif
break
;
case
IO_ERADSTICKY
:
/* sticky radiated energy */
#ifdef OUTPUTERADSTICKY
/* obsolete */
#endif
break
;
case
IO_ERADFEEDBACK
:
/* feedback injected energy */
#ifdef OUTPUTERADFEEDBACK
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
SphP
[
pindex
].
EgySpecFeedback
;
n
++
;
}
#endif
break
;
case
IO_ENERGYFLUX
:
/* energyflux */
#ifdef OUTPUTENERGYFLUX
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
SphP
[
pindex
].
EnergyFlux
;
n
++
;
}
#endif
break
;
case
IO_OPTVAR1
:
/* optional variable 1 */
#ifdef OUTPUTOPTVAR1
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
SphP
[
pindex
].
OptVar1
;
n
++
;
}
#endif
break
;
case
IO_OPTVAR2
:
/* optional variable 2 */
#ifdef OUTPUTOPTVAR2
for
(
n
=
0
;
n
<
pc
;
pindex
++
)
if
(
P
[
pindex
].
Type
==
type
)
{
*
fp
++
=
SphP
[
pindex
].
OptVar2
;
n
++
;
}
#endif
break
;
}
*
startindex
=
pindex
;
}
/*! This function tells the size of one data entry in each of the blocks
* defined for the output file. If one wants to add a new output-block, this
* function should be augmented accordingly.
*/
int
get_bytes_per_blockelement
(
enum
iofields
blocknr
)
{
int
bytes_per_blockelement
=
0
;
switch
(
blocknr
)
{
case
IO_POS
:
case
IO_VEL
:
case
IO_ACCEL
:
bytes_per_blockelement
=
3
*
sizeof
(
float
);
break
;
case
IO_ID
:
#ifdef LONGIDS
bytes_per_blockelement
=
sizeof
(
long
long
);
#else
bytes_per_blockelement
=
sizeof
(
int
);
#endif
break
;
case
IO_MASS
:
case
IO_U
:
case
IO_RHO
:
case
IO_HSML
:
case
IO_POT
:
case
IO_DTENTR
:
case
IO_TSTP
:
case
IO_ERADSPH
:
case
IO_ERADSTICKY
:
case
IO_ERADFEEDBACK
:
case
IO_ENERGYFLUX
:
case
IO_OPTVAR1
:
case
IO_OPTVAR2
:
bytes_per_blockelement
=
sizeof
(
float
);
break
;
case
IO_STAR_FORMATIONTIME
:
case
IO_INITIAL_MASS
:
case
IO_STAR_RHO
:
case
IO_STAR_HSML
:
bytes_per_blockelement
=
sizeof
(
float
);
break
;
case
IO_METALS
:
#ifdef CHIMIE
case
IO_STAR_METALS
:
bytes_per_blockelement
=
NELEMENTS
*
sizeof
(
float
);
#endif
break
;
case
IO_STAR_IDPROJ
:
#ifdef LONGIDS
bytes_per_blockelement
=
sizeof
(
long
long
);
#else
bytes_per_blockelement
=
sizeof
(
int
);
#endif
break
;
}
return
bytes_per_blockelement
;
}
/*! This function returns the type of the data contained in a given block of
* the output file. If one wants to add a new output-block, this function
* should be augmented accordingly.
*/
int
get_datatype_in_block
(
enum
iofields
blocknr
)
{
int
typekey
;
switch
(
blocknr
)
{
case
IO_ID
:
case
IO_STAR_IDPROJ
:
#ifdef LONGIDS
typekey
=
2
;
/* native long long */
#else
typekey
=
0
;
/* native int */
#endif
break
;
default:
typekey
=
1
;
/* native float */
break
;
}
return
typekey
;
}
/*! This function informs about the number of elements stored per particle for
* the given block of the output file. If one wants to add a new
* output-block, this function should be augmented accordingly.
*/
int
get_values_per_blockelement
(
enum
iofields
blocknr
)
{
int
values
=
0
;
switch
(
blocknr
)
{
case
IO_POS
:
case
IO_VEL
:
case
IO_ACCEL
:
values
=
3
;
break
;
case
IO_ID
:
case
IO_MASS
:
case
IO_U
:
case
IO_RHO
:
case
IO_HSML
:
case
IO_POT
:
case
IO_DTENTR
:
case
IO_TSTP
:
case
IO_ERADSPH
:
case
IO_ERADSTICKY
:
case
IO_ERADFEEDBACK
:
case
IO_ENERGYFLUX
:
case
IO_OPTVAR1
:
case
IO_OPTVAR2
:
case
IO_STAR_FORMATIONTIME
:
case
IO_INITIAL_MASS
:
case
IO_STAR_IDPROJ
:
case
IO_STAR_RHO
:
case
IO_STAR_HSML
:
values
=
1
;
break
;
case
IO_METALS
:
#ifdef CHIMIE
case
IO_STAR_METALS
:
values
=
NELEMENTS
;
#endif
break
;
}
return
values
;
}
/*! This function determines how many particles there are in a given block,
* based on the information in the header-structure. It also flags particle
* types that are present in the block in the typelist array. If one wants to
* add a new output-block, this function should be augmented accordingly.
*/
int
get_particles_in_block
(
enum
iofields
blocknr
,
int
*
typelist
)
{
int
i
,
nall
,
ntot_withmasses
,
ngas
,
nstars
;
nall
=
0
;
ntot_withmasses
=
0
;
for
(
i
=
0
;
i
<
6
;
i
++
)
{
typelist
[
i
]
=
0
;
if
(
header
.
npart
[
i
]
>
0
)
{
nall
+=
header
.
npart
[
i
];
typelist
[
i
]
=
1
;
}
if
(
All
.
MassTable
[
i
]
==
0
)
ntot_withmasses
+=
header
.
npart
[
i
];
}
ngas
=
header
.
npart
[
0
];
#if defined(SFR) || defined(STELLAR_PROP)
nstars
=
header
.
npart
[
ST
];
#else
nstars
=
header
.
npart
[
4
];
#endif
switch
(
blocknr
)
{
case
IO_POS
:
case
IO_VEL
:
case
IO_ACCEL
:
case
IO_TSTP
:
case
IO_ID
:
case
IO_POT
:
return
nall
;
break
;
case
IO_MASS
:
for
(
i
=
0
;
i
<
6
;
i
++
)
{
typelist
[
i
]
=
0
;
if
(
All
.
MassTable
[
i
]
==
0
&&
header
.
npart
[
i
]
>
0
)
typelist
[
i
]
=
1
;
}
return
ntot_withmasses
;
break
;
case
IO_U
:
case
IO_RHO
:
case
IO_HSML
:
case
IO_DTENTR
:
case
IO_ERADSPH
:
case
IO_ERADSTICKY
:
case
IO_ERADFEEDBACK
:
case
IO_ENERGYFLUX
:
case
IO_OPTVAR1
:
case
IO_OPTVAR2
:
case
IO_METALS
:
for
(
i
=
1
;
i
<
6
;
i
++
)
typelist
[
i
]
=
0
;
return
ngas
;
break
;
case
IO_STAR_FORMATIONTIME
:
case
IO_INITIAL_MASS
:
case
IO_STAR_IDPROJ
:
case
IO_STAR_RHO
:
case
IO_STAR_HSML
:
#if defined(STELLAR_PROP)
case
IO_STAR_METALS
:
for
(
i
=
0
;
i
<
6
;
i
++
)
{
typelist
[
i
]
=
0
;
if
(
i
==
ST
)
typelist
[
i
]
=
1
;
}
return
nstars
;
#endif
break
;
}
endrun
(
212
);
return
0
;
}
/*! This function tells whether or not a given block in the output file is
* present, depending on the type of simulation run and the compile-time
* options. If one wants to add a new output-block, this function should be
* augmented accordingly.
*/
int
blockpresent
(
enum
iofields
blocknr
)
{
#ifndef OUTPUTPOTENTIAL
if
(
blocknr
==
IO_POT
)
return
0
;
#endif
#ifndef OUTPUTACCELERATION
if
(
blocknr
==
IO_ACCEL
)
return
0
;
#endif
#ifndef OUTPUTCHANGEOFENTROPY
if
(
blocknr
==
IO_DTENTR
)
return
0
;
#endif
#ifndef OUTPUTTIMESTEP
if
(
blocknr
==
IO_TSTP
)
return
0
;
#endif
#ifndef OUTPUTERADSPH
if
(
blocknr
==
IO_ERADSPH
)
return
0
;
#endif
#ifndef OUTPUTERADSTICKY
if
(
blocknr
==
IO_ERADSTICKY
)
return
0
;
#endif
#ifndef OUTPUTERADFEEDBACK
if
(
blocknr
==
IO_ERADFEEDBACK
)
return
0
;
#endif
#ifndef OUTPUTENERGYFLUX
if
(
blocknr
==
IO_ENERGYFLUX
)
return
0
;
#endif
#ifndef OUTPUTOPTVAR1
if
(
blocknr
==
IO_OPTVAR1
)
return
0
;
#endif
#ifndef OUTPUTOPTVAR2
if
(
blocknr
==
IO_OPTVAR2
)
return
0
;
#endif
#ifndef OUTPUTSTELLAR_PROP
if
(
blocknr
==
IO_STAR_FORMATIONTIME
)
return
0
;
#endif
#ifndef OUTPUTSTELLAR_PROP
if
(
blocknr
==
IO_INITIAL_MASS
)
return
0
;
#endif
#ifndef OUTPUTSTELLAR_PROP
if
(
blocknr
==
IO_STAR_IDPROJ
)
return
0
;
#endif
#ifndef OUTPUTSTELLAR_PROP
if
(
blocknr
==
IO_STAR_RHO
)
return
0
;
#endif
#ifndef OUTPUTSTELLAR_PROP
if
(
blocknr
==
IO_STAR_HSML
)
return
0
;
#endif
#ifndef OUTPUTSTELLAR_PROP
if
(
blocknr
==
IO_STAR_METALS
)
return
0
;
#endif
#ifndef OUTPUTSTELLAR_PROP
if
(
blocknr
==
IO_METALS
)
return
0
;
#endif
return
1
;
/* default: present */
}
#ifdef BLOCK_SKIPPING
/*! This function tells whether or not a given block in the input file is
* present, depending on the type of simulation run and the compile-time
* options. If one wants to add a new input-block, this function should be
* augmented accordingly.
*/
int
blockabsent
(
enum
iofields
blocknr
)
{
/* here we will for exampe read the gas initial metallicity if needed */
#ifdef CHIMIE
#ifdef CHIMIE_INPUT_ALL
/* here, we restart from a file that contains all block*/
if
(
RestartFlag
==
0
&&
header
.
flag_metals
)
if
(
blocknr
>
IO_STAR_METALS
)
/* read all blocks */
return
1
;
else
return
0
;
#else
/* here, we restart from a file that contains only the gas metals */
if
(
RestartFlag
==
0
&&
header
.
flag_metals
)
if
(
blocknr
>
IO_METALS
)
/* read in addition IO_RHO,IO_HSML and IO_METALS */
return
1
;
else
return
0
;
#endif
#endif
if
(
RestartFlag
==
0
&&
blocknr
>
IO_U
)
return
1
;
/* ignore all other blocks in initial conditions */
return
0
;
/* default: absent */
}
#endif
/*! This function associates a short 4-character block name with each block
* number. This is stored in front of each block for snapshot
* FileFormat=2. If one wants to add a new output-block, this function should
* be augmented accordingly.
*/
void
fill_Tab_IO_Labels
(
void
)
{
enum
iofields
i
;
for
(
i
=
0
;
i
<
IO_NBLOCKS
;
i
++
)
switch
(
i
)
{
case
IO_POS
:
strncpy
(
Tab_IO_Labels
[
IO_POS
],
"POS "
,
4
);
break
;
case
IO_VEL
:
strncpy
(
Tab_IO_Labels
[
IO_VEL
],
"VEL "
,
4
);
break
;
case
IO_ID
:
strncpy
(
Tab_IO_Labels
[
IO_ID
],
"ID "
,
4
);
break
;
case
IO_MASS
:
strncpy
(
Tab_IO_Labels
[
IO_MASS
],
"MASS"
,
4
);
break
;
case
IO_U
:
strncpy
(
Tab_IO_Labels
[
IO_U
],
"U "
,
4
);
break
;
case
IO_RHO
:
strncpy
(
Tab_IO_Labels
[
IO_RHO
],
"RHO "
,
4
);
break
;
case
IO_HSML
:
strncpy
(
Tab_IO_Labels
[
IO_HSML
],
"HSML"
,
4
);
break
;
case
IO_POT
:
strncpy
(
Tab_IO_Labels
[
IO_POT
],
"POT "
,
4
);
break
;
case
IO_ACCEL
:
strncpy
(
Tab_IO_Labels
[
IO_ACCEL
],
"ACCE"
,
4
);
break
;
case
IO_DTENTR
:
strncpy
(
Tab_IO_Labels
[
IO_DTENTR
],
"ENDT"
,
4
);
break
;
case
IO_TSTP
:
strncpy
(
Tab_IO_Labels
[
IO_TSTP
],
"TSTP"
,
4
);
break
;
case
IO_ERADSPH
:
strncpy
(
Tab_IO_Labels
[
IO_ERADSPH
],
"ERADSPH"
,
4
);
break
;
case
IO_ERADSTICKY
:
strncpy
(
Tab_IO_Labels
[
IO_ERADSTICKY
],
"ERADSTICKY"
,
4
);
break
;
case
IO_ERADFEEDBACK
:
strncpy
(
Tab_IO_Labels
[
IO_ERADFEEDBACK
],
"ERADFEEDBACK"
,
4
);
break
;
case
IO_ENERGYFLUX
:
strncpy
(
Tab_IO_Labels
[
IO_ENERGYFLUX
],
"ENERGYFLUX"
,
4
);
break
;
case
IO_OPTVAR1
:
strncpy
(
Tab_IO_Labels
[
IO_OPTVAR1
],
"OPTVAR1"
,
4
);
break
;
case
IO_OPTVAR2
:
strncpy
(
Tab_IO_Labels
[
IO_OPTVAR2
],
"OPTVAR2"
,
4
);
break
;
case
IO_STAR_FORMATIONTIME
:
strncpy
(
Tab_IO_Labels
[
IO_STAR_FORMATIONTIME
],
"STAR_FORMATIONTIME"
,
4
);
break
;
case
IO_INITIAL_MASS
:
strncpy
(
Tab_IO_Labels
[
IO_INITIAL_MASS
],
"INITIAL_MASS"
,
4
);
break
;
case
IO_STAR_IDPROJ
:
strncpy
(
Tab_IO_Labels
[
IO_STAR_IDPROJ
],
"STAR_IDPROJ"
,
4
);
break
;
case
IO_STAR_RHO
:
strncpy
(
Tab_IO_Labels
[
IO_STAR_RHO
],
"STAR_RHO"
,
4
);
break
;
case
IO_STAR_HSML
:
strncpy
(
Tab_IO_Labels
[
IO_STAR_HSML
],
"STAR_HSML"
,
4
);
break
;
case
IO_STAR_METALS
:
strncpy
(
Tab_IO_Labels
[
IO_STAR_METALS
],
"STAR_METALS"
,
4
);
break
;
case
IO_METALS
:
strncpy
(
Tab_IO_Labels
[
IO_METALS
],
"METALS"
,
4
);
break
;
}
}
/*! This function returns a descriptive character string that describes the
* name of the block when the HDF5 file format is used. If one wants to add
* a new output-block, this function should be augmented accordingly.
*/
void
get_dataset_name
(
enum
iofields
blocknr
,
char
*
buf
)
{
strcpy
(
buf
,
"default"
);
switch
(
blocknr
)
{
case
IO_POS
:
strcpy
(
buf
,
"Coordinates"
);
break
;
case
IO_VEL
:
strcpy
(
buf
,
"Velocities"
);
break
;
case
IO_ID
:
strcpy
(
buf
,
"ParticleIDs"
);
break
;
case
IO_MASS
:
strcpy
(
buf
,
"Masses"
);
break
;
case
IO_U
:
strcpy
(
buf
,
"InternalEnergy"
);
break
;
case
IO_RHO
:
strcpy
(
buf
,
"Density"
);
break
;
case
IO_HSML
:
strcpy
(
buf
,
"SmoothingLength"
);
break
;
case
IO_POT
:
strcpy
(
buf
,
"Potential"
);
break
;
case
IO_ACCEL
:
strcpy
(
buf
,
"Acceleration"
);
break
;
case
IO_DTENTR
:
strcpy
(
buf
,
"RateOfChangeOfEntropy"
);
break
;
case
IO_TSTP
:
strcpy
(
buf
,
"TimeStep"
);
break
;
case
IO_ERADSPH
:
strcpy
(
buf
,
"EnergyRadiatedSph"
);
break
;
case
IO_ERADSTICKY
:
strcpy
(
buf
,
"EnergyRadiatedSticky"
);
break
;
case
IO_ERADFEEDBACK
:
strcpy
(
buf
,
"EnergyRadiatedFeedback"
);
break
;
case
IO_ENERGYFLUX
:
strcpy
(
buf
,
"EnergyFlux"
);
break
;
case
IO_OPTVAR1
:
strcpy
(
buf
,
"OptVar1"
);
break
;
case
IO_OPTVAR2
:
strcpy
(
buf
,
"OptVar2"
);
break
;
case
IO_STAR_FORMATIONTIME
:
strcpy
(
buf
,
"StarFormationTime"
);
break
;
case
IO_INITIAL_MASS
:
strcpy
(
buf
,
"InitialMass"
);
break
;
case
IO_STAR_IDPROJ
:
strcpy
(
buf
,
"StarIDProj"
);
break
;
case
IO_STAR_RHO
:
strcpy
(
buf
,
"StarRho"
);
break
;
case
IO_STAR_HSML
:
strcpy
(
buf
,
"StarHsml"
);
break
;
case
IO_STAR_METALS
:
strcpy
(
buf
,
"StarMetals"
);
break
;
case
IO_METALS
:
strcpy
(
buf
,
"Metals"
);
break
;
}
}
/*! This function writes an actual snapshot file containing the data from
* processors 'writeTask' to 'lastTask'. 'writeTask' is the one that actually
* writes. Each snapshot file contains a header first, then particle
* positions, velocities and ID's. Particle masses are written only for
* those particle types with zero entry in MassTable. After that, first the
* internal energies u, and then the density is written for the SPH
* particles. If cooling is enabled, mean molecular weight and neutral
* hydrogen abundance are written for the gas particles. This is followed by
* the SPH smoothing length and further blocks of information, depending on
* included physics and compile-time flags. If HDF5 is used, the header is
* stored in a group called "/Header", and the particle data is stored
* separately for each particle type in groups calles "/PartType0",
* "/PartType1", etc. The sequence of the blocks is unimportant in this case.
*/
void
write_file
(
char
*
fname
,
int
writeTask
,
int
lastTask
)
{
int
type
,
bytes_per_blockelement
,
npart
,
nextblock
,
typelist
[
6
];
int
n_for_this_task
,
ntask
,
n
,
p
,
pc
,
offset
=
0
,
task
;
int
blockmaxlen
,
ntot_type
[
6
],
nn
[
6
];
enum
iofields
blocknr
;
int
blksize
;
int
i
;
MPI_Status
status
;
FILE
*
fd
=
0
;
#ifdef HAVE_HDF5
hid_t
hdf5_file
=
0
,
hdf5_grp
[
6
],
hdf5_headergrp
=
0
,
hdf5_dataspace_memory
;
hid_t
hdf5_datatype
=
0
,
hdf5_dataspace_in_file
=
0
,
hdf5_dataset
=
0
;
herr_t
hdf5_status
;
hsize_t
dims
[
2
],
count
[
2
],
start
[
2
];
int
rank
,
pcsum
=
0
;
char
buf
[
500
];
#endif
#define SKIP {my_fwrite(&blksize,sizeof(int),1,fd);}
/* determine particle numbers of each type in file */
if
(
ThisTask
==
writeTask
)
{
for
(
n
=
0
;
n
<
6
;
n
++
)
ntot_type
[
n
]
=
n_type
[
n
];
for
(
task
=
writeTask
+
1
;
task
<=
lastTask
;
task
++
)
{
MPI_Recv
(
&
nn
[
0
],
6
,
MPI_INT
,
task
,
TAG_LOCALN
,
MPI_COMM_WORLD
,
&
status
);
for
(
n
=
0
;
n
<
6
;
n
++
)
ntot_type
[
n
]
+=
nn
[
n
];
}
for
(
task
=
writeTask
+
1
;
task
<=
lastTask
;
task
++
)
MPI_Send
(
&
ntot_type
[
0
],
6
,
MPI_INT
,
task
,
TAG_N
,
MPI_COMM_WORLD
);
}
else
{
MPI_Send
(
&
n_type
[
0
],
6
,
MPI_INT
,
writeTask
,
TAG_LOCALN
,
MPI_COMM_WORLD
);
MPI_Recv
(
&
ntot_type
[
0
],
6
,
MPI_INT
,
writeTask
,
TAG_N
,
MPI_COMM_WORLD
,
&
status
);
}
/* fill file header */
for
(
n
=
0
;
n
<
6
;
n
++
)
{
header
.
npart
[
n
]
=
ntot_type
[
n
];
header
.
npartTotal
[
n
]
=
(
unsigned
int
)
ntot_type_all
[
n
];
header
.
npartTotalHighWord
[
n
]
=
(
unsigned
int
)
(
ntot_type_all
[
n
]
>>
32
);
}
for
(
n
=
0
;
n
<
6
;
n
++
)
header
.
mass
[
n
]
=
All
.
MassTable
[
n
];
header
.
time
=
All
.
Time
;
if
(
All
.
ComovingIntegrationOn
)
header
.
redshift
=
1.0
/
All
.
Time
-
1
;
else
header
.
redshift
=
0
;
header
.
flag_sfr
=
0
;
header
.
flag_feedback
=
0
;
header
.
flag_cooling
=
0
;
header
.
flag_stellarage
=
0
;
header
.
flag_metals
=
0
;
#ifdef MULTIPHASE
header
.
critical_energy_spec
=
All
.
CriticalEgySpec
;
#endif
#ifdef COOLING
header
.
flag_cooling
=
1
;
#endif
#ifdef SFR
header
.
flag_sfr
=
1
;
#ifdef OUTPUTSTELLAR_PROP
header
.
flag_stellarage
=
1
;
header
.
flag_metals
=
NELEMENTS
;
#endif
#ifdef FEEDBACK
header
.
flag_feedback
=
1
;
#endif
#endif
header
.
num_files
=
All
.
NumFilesPerSnapshot
;
header
.
BoxSize
=
All
.
BoxSize
;
header
.
Omega0
=
All
.
Omega0
;
header
.
OmegaLambda
=
All
.
OmegaLambda
;
header
.
HubbleParam
=
All
.
HubbleParam
;
/* set fill to " " : yr Thu Aug 13 17:34:07 CEST 2009*/
memset
(
header
.
fill
,
' '
,
sizeof
(
header
.
fill
));
/* fill file chimie-header */
header
.
flag_chimie_extraheader
=
0
;
#ifdef CHIMIE_EXTRAHEADER
header
.
flag_chimie_extraheader
=
1
;
chimie_extraheader
.
nelts
=
get_nelts
();
for
(
i
=
0
;
i
<
get_nelts
();
i
++
)
{
chimie_extraheader
.
SolarAbundances
[
i
]
=
get_SolarAbundance
(
i
);
}
memset
(
chimie_extraheader
.
labels
,
' '
,
sizeof
(
chimie_extraheader
.
labels
));
for
(
i
=
0
,
n
=
0
;
i
<
get_nelts
();
i
++
)
{
strcpy
(
&
chimie_extraheader
.
labels
[
n
],
get_Element
(
i
));
n
+=
strlen
(
get_Element
(
i
));
strncpy
(
&
chimie_extraheader
.
labels
[
n
++
],
","
,
1
);
}
#endif
/* open file and write header */
if
(
ThisTask
==
writeTask
)
{
if
(
All
.
SnapFormat
==
3
)
{
#ifdef HAVE_HDF5
sprintf
(
buf
,
"%s.hdf5"
,
fname
);
hdf5_file
=
H5Fcreate
(
buf
,
H5F_ACC_TRUNC
,
H5P_DEFAULT
,
H5P_DEFAULT
);
hdf5_headergrp
=
H5Gcreate
(
hdf5_file
,
"/Header"
,
0
);
for
(
type
=
0
;
type
<
6
;
type
++
)
{
if
(
header
.
npart
[
type
]
>
0
)
{
sprintf
(
buf
,
"/PartType%d"
,
type
);
hdf5_grp
[
type
]
=
H5Gcreate
(
hdf5_file
,
buf
,
0
);
}
}
write_header_attributes_in_hdf5
(
hdf5_headergrp
);
#endif
}
else
{
if
(
!
(
fd
=
fopen
(
fname
,
"w"
)))
{
printf
(
"can't open file `%s' for writing snapshot.
\n
"
,
fname
);
endrun
(
123
);
}
if
(
All
.
SnapFormat
==
2
)
{
blksize
=
sizeof
(
int
)
+
4
*
sizeof
(
char
);
SKIP
;
my_fwrite
(
"HEAD"
,
sizeof
(
char
),
4
,
fd
);
nextblock
=
sizeof
(
header
)
+
2
*
sizeof
(
int
);
my_fwrite
(
&
nextblock
,
sizeof
(
int
),
1
,
fd
);
SKIP
;
}
blksize
=
sizeof
(
header
);
SKIP
;
my_fwrite
(
&
header
,
sizeof
(
header
),
1
,
fd
);
SKIP
;
#ifdef CHIMIE_EXTRAHEADER
blksize
=
sizeof
(
chimie_extraheader
);
SKIP
;
my_fwrite
(
&
chimie_extraheader
,
sizeof
(
chimie_extraheader
),
1
,
fd
);
SKIP
;
#endif
}
}
ntask
=
lastTask
-
writeTask
+
1
;
for
(
blocknr
=
0
;
blocknr
<
IO_NBLOCKS
;
blocknr
++
)
{
if
(
blockpresent
(
blocknr
))
{
bytes_per_blockelement
=
get_bytes_per_blockelement
(
blocknr
);
blockmaxlen
=
((
int
)
(
All
.
BufferSize
*
1024
*
1024
))
/
bytes_per_blockelement
;
npart
=
get_particles_in_block
(
blocknr
,
&
typelist
[
0
]);
if
(
npart
>
0
)
{
if
(
ThisTask
==
writeTask
)
{
if
(
All
.
SnapFormat
==
1
||
All
.
SnapFormat
==
2
)
{
if
(
All
.
SnapFormat
==
2
)
{
blksize
=
sizeof
(
int
)
+
4
*
sizeof
(
char
);
SKIP
;
my_fwrite
(
Tab_IO_Labels
[
blocknr
],
sizeof
(
char
),
4
,
fd
);
nextblock
=
npart
*
bytes_per_blockelement
+
2
*
sizeof
(
int
);
my_fwrite
(
&
nextblock
,
sizeof
(
int
),
1
,
fd
);
SKIP
;
}
blksize
=
npart
*
bytes_per_blockelement
;
SKIP
;
}
}
for
(
type
=
0
;
type
<
6
;
type
++
)
{
if
(
typelist
[
type
])
{
#ifdef HAVE_HDF5
if
(
ThisTask
==
writeTask
&&
All
.
SnapFormat
==
3
&&
header
.
npart
[
type
]
>
0
)
{
switch
(
get_datatype_in_block
(
blocknr
))
{
case
0
:
hdf5_datatype
=
H5Tcopy
(
H5T_NATIVE_UINT
);
break
;
case
1
:
hdf5_datatype
=
H5Tcopy
(
H5T_NATIVE_FLOAT
);
break
;
case
2
:
hdf5_datatype
=
H5Tcopy
(
H5T_NATIVE_UINT64
);
break
;
}
dims
[
0
]
=
header
.
npart
[
type
];
dims
[
1
]
=
get_values_per_blockelement
(
blocknr
);
if
(
dims
[
1
]
==
1
)
rank
=
1
;
else
rank
=
2
;
get_dataset_name
(
blocknr
,
buf
);
hdf5_dataspace_in_file
=
H5Screate_simple
(
rank
,
dims
,
NULL
);
hdf5_dataset
=
H5Dcreate
(
hdf5_grp
[
type
],
buf
,
hdf5_datatype
,
hdf5_dataspace_in_file
,
H5P_DEFAULT
);
pcsum
=
0
;
}
#endif
for
(
task
=
writeTask
,
offset
=
0
;
task
<=
lastTask
;
task
++
)
{
if
(
task
==
ThisTask
)
{
n_for_this_task
=
n_type
[
type
];
for
(
p
=
writeTask
;
p
<=
lastTask
;
p
++
)
if
(
p
!=
ThisTask
)
MPI_Send
(
&
n_for_this_task
,
1
,
MPI_INT
,
p
,
TAG_NFORTHISTASK
,
MPI_COMM_WORLD
);
}
else
MPI_Recv
(
&
n_for_this_task
,
1
,
MPI_INT
,
task
,
TAG_NFORTHISTASK
,
MPI_COMM_WORLD
,
&
status
);
while
(
n_for_this_task
>
0
)
{
pc
=
n_for_this_task
;
if
(
pc
>
blockmaxlen
)
pc
=
blockmaxlen
;
if
(
ThisTask
==
task
)
fill_write_buffer
(
blocknr
,
&
offset
,
pc
,
type
);
if
(
ThisTask
==
writeTask
&&
task
!=
writeTask
)
MPI_Recv
(
CommBuffer
,
bytes_per_blockelement
*
pc
,
MPI_BYTE
,
task
,
TAG_PDATA
,
MPI_COMM_WORLD
,
&
status
);
if
(
ThisTask
!=
writeTask
&&
task
==
ThisTask
)
MPI_Ssend
(
CommBuffer
,
bytes_per_blockelement
*
pc
,
MPI_BYTE
,
writeTask
,
TAG_PDATA
,
MPI_COMM_WORLD
);
if
(
ThisTask
==
writeTask
)
{
if
(
All
.
SnapFormat
==
3
)
{
#ifdef HAVE_HDF5
start
[
0
]
=
pcsum
;
start
[
1
]
=
0
;
count
[
0
]
=
pc
;
count
[
1
]
=
get_values_per_blockelement
(
blocknr
);
pcsum
+=
pc
;
H5Sselect_hyperslab
(
hdf5_dataspace_in_file
,
H5S_SELECT_SET
,
start
,
NULL
,
count
,
NULL
);
dims
[
0
]
=
pc
;
dims
[
1
]
=
get_values_per_blockelement
(
blocknr
);
hdf5_dataspace_memory
=
H5Screate_simple
(
rank
,
dims
,
NULL
);
hdf5_status
=
H5Dwrite
(
hdf5_dataset
,
hdf5_datatype
,
hdf5_dataspace_memory
,
hdf5_dataspace_in_file
,
H5P_DEFAULT
,
CommBuffer
);
H5Sclose
(
hdf5_dataspace_memory
);
#endif
}
else
my_fwrite
(
CommBuffer
,
bytes_per_blockelement
,
pc
,
fd
);
}
n_for_this_task
-=
pc
;
}
}
#ifdef HAVE_HDF5
if
(
ThisTask
==
writeTask
&&
All
.
SnapFormat
==
3
&&
header
.
npart
[
type
]
>
0
)
{
if
(
All
.
SnapFormat
==
3
)
{
H5Dclose
(
hdf5_dataset
);
H5Sclose
(
hdf5_dataspace_in_file
);
H5Tclose
(
hdf5_datatype
);
}
}
#endif
}
}
if
(
ThisTask
==
writeTask
)
{
if
(
All
.
SnapFormat
==
1
||
All
.
SnapFormat
==
2
)
SKIP
;
}
}
}
}
if
(
ThisTask
==
writeTask
)
{
if
(
All
.
SnapFormat
==
3
)
{
#ifdef HAVE_HDF5
for
(
type
=
5
;
type
>=
0
;
type
--
)
if
(
header
.
npart
[
type
]
>
0
)
H5Gclose
(
hdf5_grp
[
type
]);
H5Gclose
(
hdf5_headergrp
);
H5Fclose
(
hdf5_file
);
#endif
}
else
fclose
(
fd
);
}
}
/*! This function writes the header information in case HDF5 is selected as
* file format.
*/
#ifdef HAVE_HDF5
void
write_header_attributes_in_hdf5
(
hid_t
handle
)
{
hsize_t
adim
[
1
]
=
{
6
};
hid_t
hdf5_dataspace
,
hdf5_attribute
;
hdf5_dataspace
=
H5Screate
(
H5S_SIMPLE
);
H5Sset_extent_simple
(
hdf5_dataspace
,
1
,
adim
,
NULL
);
hdf5_attribute
=
H5Acreate
(
handle
,
"NumPart_ThisFile"
,
H5T_NATIVE_INT
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_UINT
,
header
.
npart
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SIMPLE
);
H5Sset_extent_simple
(
hdf5_dataspace
,
1
,
adim
,
NULL
);
hdf5_attribute
=
H5Acreate
(
handle
,
"NumPart_Total"
,
H5T_NATIVE_UINT
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_UINT
,
header
.
npartTotal
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SIMPLE
);
H5Sset_extent_simple
(
hdf5_dataspace
,
1
,
adim
,
NULL
);
hdf5_attribute
=
H5Acreate
(
handle
,
"NumPart_Total_HW"
,
H5T_NATIVE_UINT
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_UINT
,
header
.
npartTotalHighWord
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SIMPLE
);
H5Sset_extent_simple
(
hdf5_dataspace
,
1
,
adim
,
NULL
);
hdf5_attribute
=
H5Acreate
(
handle
,
"MassTable"
,
H5T_NATIVE_DOUBLE
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_DOUBLE
,
header
.
mass
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SCALAR
);
hdf5_attribute
=
H5Acreate
(
handle
,
"Time"
,
H5T_NATIVE_DOUBLE
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_DOUBLE
,
&
header
.
time
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SCALAR
);
hdf5_attribute
=
H5Acreate
(
handle
,
"Redshift"
,
H5T_NATIVE_DOUBLE
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_DOUBLE
,
&
header
.
redshift
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SCALAR
);
hdf5_attribute
=
H5Acreate
(
handle
,
"BoxSize"
,
H5T_NATIVE_DOUBLE
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_DOUBLE
,
&
header
.
BoxSize
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SCALAR
);
hdf5_attribute
=
H5Acreate
(
handle
,
"NumFilesPerSnapshot"
,
H5T_NATIVE_INT
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_INT
,
&
header
.
num_files
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SCALAR
);
hdf5_attribute
=
H5Acreate
(
handle
,
"Omega0"
,
H5T_NATIVE_DOUBLE
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_DOUBLE
,
&
header
.
Omega0
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SCALAR
);
hdf5_attribute
=
H5Acreate
(
handle
,
"OmegaLambda"
,
H5T_NATIVE_DOUBLE
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_DOUBLE
,
&
header
.
OmegaLambda
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SCALAR
);
hdf5_attribute
=
H5Acreate
(
handle
,
"HubbleParam"
,
H5T_NATIVE_DOUBLE
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_DOUBLE
,
&
header
.
HubbleParam
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SCALAR
);
hdf5_attribute
=
H5Acreate
(
handle
,
"Flag_Sfr"
,
H5T_NATIVE_INT
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_INT
,
&
header
.
flag_sfr
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SCALAR
);
hdf5_attribute
=
H5Acreate
(
handle
,
"Flag_Cooling"
,
H5T_NATIVE_INT
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_INT
,
&
header
.
flag_cooling
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SCALAR
);
hdf5_attribute
=
H5Acreate
(
handle
,
"Flag_StellarAge"
,
H5T_NATIVE_INT
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_INT
,
&
header
.
flag_stellarage
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SCALAR
);
hdf5_attribute
=
H5Acreate
(
handle
,
"Flag_Metals"
,
H5T_NATIVE_INT
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_INT
,
&
header
.
flag_metals
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
hdf5_dataspace
=
H5Screate
(
H5S_SCALAR
);
hdf5_attribute
=
H5Acreate
(
handle
,
"Flag_Feedback"
,
H5T_NATIVE_INT
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_INT
,
&
header
.
flag_feedback
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
header
.
flag_entropy_instead_u
=
0
;
hdf5_dataspace
=
H5Screate
(
H5S_SIMPLE
);
H5Sset_extent_simple
(
hdf5_dataspace
,
1
,
adim
,
NULL
);
hdf5_attribute
=
H5Acreate
(
handle
,
"Flag_Entropy_ICs"
,
H5T_NATIVE_UINT
,
hdf5_dataspace
,
H5P_DEFAULT
);
H5Awrite
(
hdf5_attribute
,
H5T_NATIVE_UINT
,
header
.
flag_entropy_instead_u
);
H5Aclose
(
hdf5_attribute
);
H5Sclose
(
hdf5_dataspace
);
}
#endif
/*! This catches I/O errors occuring for my_fwrite(). In this case we
* better stop.
*/
size_t
my_fwrite
(
void
*
ptr
,
size_t
size
,
size_t
nmemb
,
FILE
*
stream
)
{
size_t
nwritten
;
if
((
nwritten
=
fwrite
(
ptr
,
size
,
nmemb
,
stream
))
!=
nmemb
)
{
printf
(
"I/O error (fwrite) on task=%d has occured: %s
\n
"
,
ThisTask
,
strerror
(
errno
));
fflush
(
stdout
);
endrun
(
777
);
}
return
nwritten
;
}
/*! This catches I/O errors occuring for fread(). In this case we
* better stop.
*/
size_t
my_fread
(
void
*
ptr
,
size_t
size
,
size_t
nmemb
,
FILE
*
stream
)
{
size_t
nread
;
if
((
nread
=
fread
(
ptr
,
size
,
nmemb
,
stream
))
!=
nmemb
)
{
printf
(
"I/O error (fread) on task=%d has occured: %s
\n
"
,
ThisTask
,
strerror
(
errno
));
fflush
(
stdout
);
endrun
(
778
);
}
return
nread
;
}
Event Timeline
Log In to Comment