aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--configuration.ccl3
-rw-r--r--interface.ccl5
-rw-r--r--param.ccl13
-rw-r--r--src/multipole.cc29
-rw-r--r--src/utils.cc120
-rw-r--r--src/utils.hh3
6 files changed, 165 insertions, 8 deletions
diff --git a/configuration.ccl b/configuration.ccl
new file mode 100644
index 0000000..b3f112e
--- /dev/null
+++ b/configuration.ccl
@@ -0,0 +1,3 @@
+OPTIONAL HDF5
+{
+}
diff --git a/interface.ccl b/interface.ccl
index 463a8ce..a8eda59 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,6 +2,11 @@
implements: multipole
inherits: Grid
+CCTK_INT FUNCTION IO_TruncateOutputFiles \
+ (CCTK_POINTER_TO_CONST IN GH)
+
+REQUIRES FUNCTION IO_TruncateOutputFiles
+
public:
CCTK_REAL harmonics type=GF timelevels=1
diff --git a/param.ccl b/param.ccl
index 2f089dc..9af3650 100644
--- a/param.ccl
+++ b/param.ccl
@@ -137,3 +137,16 @@ CCTK_INT m_mode "all modes: Up to which m mode to calculate/ specific mode: whic
{
-100:* :: "Deprecated"
} -100
+
+CCTK_BOOLEAN output_ascii "Output a simple ASCII file for each mode at each radius"
+{
+} "yes"
+
+CCTK_BOOLEAN output_hdf5 "Output an HDF5 file for each variable containing one dataset per mode at each radius"
+{
+} "no"
+
+CCTK_INT hdf5_chunk_size "How many iterations to preallocate in extensible HDF5 datasets"
+{
+ 1:* :: "Any integer"
+} 200
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[],
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-10 21:46:28 +0000