diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/multipole.cc | 29 | ||||
| -rw-r--r-- | src/utils.cc | 120 | ||||
| -rw-r--r-- | src/utils.hh | 3 |
3 files changed, 144 insertions, 8 deletions
diff --git a/src/multipole.cc b/src/multipole.cc index 975d589..8a6a6ea 100644 --- a/src/multipole.cc +++ b/src/multipole.cc @@ -97,15 +97,28 @@ static void parse_variables_string(const string &var_string, variable_desc v[max static void output_mode(CCTK_ARGUMENTS, const variable_desc *v, CCTK_REAL rad, CCTK_REAL real_lm, CCTK_REAL imag_lm, int l, int m) { - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; - if (CCTK_MyProc(cctkGH) == 0) + if (output_ascii) + { + if (CCTK_MyProc(cctkGH) == 0) + { + ostringstream name; + name << "mp_" << v->name << "_l" << l << "_m" << m << + "_r" << setiosflags(ios::fixed) << setprecision(2) << rad << ".asc"; + Multipole_OutputComplexToFile(CCTK_PASS_CTOC, name.str(), real_lm, imag_lm); + } + } + if (output_hdf5) { - ostringstream name; - name << "mp_" << v->name << "_l" << l << "_m" << m << - "_r" << setiosflags(ios::fixed) << setprecision(2) << rad << ".asc"; - Multipole_OutputComplexToFile(CCTK_PASS_CTOC, name.str(), real_lm, imag_lm); + if (CCTK_MyProc(cctkGH) == 0) + { + ostringstream h5datasetname; + h5datasetname << "l" << l << "_m" << m << "_r" << setiosflags(ios::fixed) << setprecision(2) << rad; + Multipole_OutputComplexToH5File(CCTK_PASS_CTOC, "mp_" + v->name + ".h5", h5datasetname.str(), + real_lm, imag_lm); + } } } @@ -116,7 +129,7 @@ static void output_1D(CCTK_ARGUMENTS, const variable_desc *v, CCTK_REAL rad, DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - if (CCTK_MyProc(cctkGH) == 0) + if (CCTK_MyProc(cctkGH) == 0 && output_ascii) { if (out_1d_every != 0 && cctk_iteration % out_1d_every == 0) { diff --git a/src/utils.cc b/src/utils.cc index 168acc2..7140718 100644 --- a/src/utils.cc +++ b/src/utils.cc @@ -2,12 +2,25 @@ #include <string.h> #include <math.h> #include <string> +#include <sys/stat.h> +#include <errno.h> #include <iostream> #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" +#ifdef HAVE_CAPABILITY_HDF5 +// We currently support the HDF5 1.6 API (and when using 1.8 the +// compatibility mode introduced by H5_USE_16_API). Several machines +// in SimFactory use HDF5 1.6, so we cannot drop support for it. It +// seems it is hard to support both the 1.6 and 1.8 API +// simultaneously; for example H5Fopen takes a different number of +// arguments in the two versions. +#define H5_USE_16_API +#include <hdf5.h> +#endif + #include "utils.hh" #include "integrate.hh" @@ -133,6 +146,113 @@ void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, CCTK_REAL } } +//////////////////////////////////////////////////////////////////////// +// HDF5 complex output +//////////////////////////////////////////////////////////////////////// + +#ifdef HAVE_CAPABILITY_HDF5 + +bool dataset_exists(hid_t file, const string &dataset_name) +{ + // To test whether a dataset exists, the recommended way in API 1.6 + // is to use H5Gget_objinfo, but this prints an error to stderr if + // the dataset does not exist. We explicitly avoid this by wrapping + // the call in H5E_BEGIN_TRY/H5E_END_TRY statements. In 1.8, + // H5Gget_objinfo is deprecated, and H5Lexists does the job. See + // http://www.mail-archive.com/hdf-forum@hdfgroup.org/msg00125.html + + #if 1 + bool exists; + H5E_BEGIN_TRY + { + exists = H5Gget_objinfo(file, dataset_name.c_str(), 1, NULL) >= 0; + } + H5E_END_TRY; + return exists; + #else + return H5Lexists(file, dataset_name.c_str(), H5P_DEFAULT); + #endif +} + +void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const string &basename, const string &datasetname, + CCTK_REAL redata, CCTK_REAL imdata) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + string output_name = out_dir + string("/") + basename; + int status = 0; + static int first_time = true; + + hid_t file = + (first_time && IO_TruncateOutputFiles(cctkGH)) ? + H5Fcreate(output_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT) : + H5Fopen(output_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT); + + first_time = false; + + hid_t dataset = -1; + + if (dataset_exists(file, datasetname)) + { + dataset = H5Dopen(file, datasetname.c_str()); + } + else + { + hsize_t dims[2] = {0,3}; + hsize_t maxdims[2] = {H5S_UNLIMITED, 3}; + hid_t dataspace = H5Screate_simple(2, dims, maxdims); + + hid_t cparms = -1; + hsize_t chunk_dims[2] = {hdf5_chunk_size,3}; + cparms = H5Pcreate (H5P_DATASET_CREATE); + status = H5Pset_chunk(cparms, 2, chunk_dims); + + dataset = H5Dcreate(file, datasetname.c_str(), H5T_NATIVE_DOUBLE, dataspace, cparms); + H5Pclose(cparms); + } + + hid_t filespace = H5Dget_space (dataset); + hsize_t filedims[2]; + hsize_t maxdims[2]; + status = H5Sget_simple_extent_dims(filespace, filedims, maxdims); + + filedims[0] += 1; + hsize_t size[2] = {filedims[0], filedims[1]}; + status = H5Dextend (dataset, size); + H5Sclose(filespace); + + /* Select a hyperslab */ + hsize_t offset[2] = {filedims[0]-1, 0}; + hsize_t dims2[2] = {1,3}; + filespace = H5Dget_space (dataset); + status = H5Sselect_hyperslab (filespace, H5S_SELECT_SET, offset, NULL, + dims2, NULL); + + CCTK_REAL data[] = {cctk_time, redata, imdata}; + + hid_t memdataspace = H5Screate_simple(2, dims2, NULL); + + /* Write the data to the hyperslab */ + status = H5Dwrite (dataset, H5T_NATIVE_DOUBLE, memdataspace, filespace, + H5P_DEFAULT, data); + + H5Dclose(dataset); + H5Sclose(filespace); + H5Sclose(memdataspace); + + status = H5Fclose(file); +} + +#else + +void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const string &basename, const string &datasetname, + CCTK_REAL redata, CCTK_REAL imdata) +{ + CCTK_WARN(0,"HDF5 output has been requested but Cactus has been compiled without HDF5 support"); +} + +#endif //////////////////////////////////////////////////////////////////////// // Coordinates diff --git a/src/utils.hh b/src/utils.hh index 51c76a7..7d10e8b 100644 --- a/src/utils.hh +++ b/src/utils.hh @@ -20,6 +20,9 @@ void Multipole_Output1D(CCTK_ARGUMENTS, const string &name, int array_size, void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, CCTK_REAL redata, CCTK_REAL imdata); +void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const string &filename, const string &dataset, + CCTK_REAL redata, CCTK_REAL imdata); + void Multipole_CoordSetup(int ntheta, int nphi, CCTK_REAL xhat[], CCTK_REAL yhat[], CCTK_REAL zhat[], CCTK_REAL th[], |