diff --git a/bfps/cpp/field.cpp b/bfps/cpp/field.cpp
index d5696d56cce92fac5d78ea80120e400677f1e38b..e64a1b95a365af3adfe3bb92aab252f83e2fa27d 100644
--- a/bfps/cpp/field.cpp
+++ b/bfps/cpp/field.cpp
@@ -612,6 +612,197 @@ int field<rnumber, be, fc>::write_0slice(
return EXIT_SUCCESS;
}
+template <typename rnumber,
+ field_backend be,
+ field_components fc>
+int field<rnumber, be, fc>::write_filtered(
+ const std::string fname,
+ const std::string field_name,
+ const int iteration,
+ int nx,
+ int ny,
+ int nz)
+{
+ /* file dataset has same dimensions as field */
+ TIMEZONE("field::write_filtered");
+ // only works in Fourier representation
+ assert(!this->real_space_representation);
+ assert(hsize_t(nx) <= this->rlayout->sizes[2]);
+ assert(hsize_t(ny) <= this->rlayout->sizes[1]);
+ assert(hsize_t(nz) <= this->rlayout->sizes[0]);
+ // current algorithm only works for more than one process
+ assert(this->nprocs >= 2);
+ hid_t file_id, dset_id, plist_id;
+ dset_id = H5I_BADID;
+ std::string dset_name = (
+ "/" + field_name +
+ "/complex" +
+ "/" + std::to_string(iteration));
+
+ /* open/create file */
+ plist_id = H5Pcreate(H5P_FILE_ACCESS);
+ H5Pset_fapl_mpio(plist_id, this->comm, MPI_INFO_NULL);
+ bool file_exists = false;
+ struct stat file_buffer;
+ file_exists = (stat(fname.c_str(), &file_buffer) == 0);
+ if (file_exists)
+ file_id = H5Fopen(fname.c_str(), H5F_ACC_RDWR, plist_id);
+ else
+ file_id = H5Fcreate(fname.c_str(), H5F_ACC_EXCL, H5P_DEFAULT, plist_id);
+ assert(file_id >= 0);
+ H5Pclose(plist_id);
+
+ /* generic space initialization */
+ hid_t fspace, mspace;
+ hsize_t count[ndim(fc)], offset[ndim(fc)], dims[ndim(fc)], fdims[ndim(fc)];
+ hsize_t memoffset[ndim(fc)], memshape[ndim(fc)];
+
+ // set up dimensions
+ for (unsigned int i=3; i<ndim(fc); i++)
+ {
+ count [i] = this->clayout->subsizes[i];
+ offset[i] = this->clayout->starts[i];
+ dims [i] = this->clayout->sizes[i];
+ memshape [i] = count[i];
+ memoffset[i] = 0;
+ }
+ // these are dimensions of dataset, needed
+ // to create dataset
+ //dims[0] = nz;
+ dims[0] = ny;
+ dims[1] = nz;
+ dims[2] = nx/2+1;
+
+ /* open/create data set */
+ if (!H5Lexists(file_id, field_name.c_str(), H5P_DEFAULT))
+ {
+ hid_t gid_tmp = H5Gcreate(
+ file_id, field_name.c_str(),
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ H5Gclose(gid_tmp);
+ }
+ if (!H5Lexists(file_id, (field_name + "/complex").c_str(), H5P_DEFAULT))
+ {
+ hid_t gid_tmp = H5Gcreate(
+ file_id, ("/" + field_name + "/complex").c_str(),
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ H5Gclose(gid_tmp);
+ }
+ if (H5Lexists(file_id, dset_name.c_str(), H5P_DEFAULT))
+ {
+ dset_id = H5Dopen(file_id, dset_name.c_str(), H5P_DEFAULT);
+ fspace = H5Dget_space(dset_id);
+ }
+ else
+ {
+ fspace = H5Screate_simple(
+ ndim(fc),
+ dims,
+ NULL);
+ /* chunking needs to go in here */
+ dset_id = H5Dcreate(
+ file_id,
+ dset_name.c_str(),
+ this->cnumber_H5T,
+ fspace,
+ H5P_DEFAULT,
+ H5P_DEFAULT,
+ H5P_DEFAULT);
+ }
+ /* check file space */
+ int ndims_fspace = H5Sget_simple_extent_dims(fspace, fdims, NULL);
+ assert(((unsigned int)(ndims_fspace)) == ndim(fc));
+ for (unsigned int i=0; i<ndim(fc); i++)
+ {
+ assert(dims[i] == fdims[i]);
+ }
+ /* both dset_id and fspace now have sane values */
+
+ /// set up counts and offsets
+ /// x is easy, since only positive modes are present
+ count [2] = nx/2+1;
+ offset[2] = 0;
+ memshape [2] = this->clayout->subsizes[2];
+ memoffset[2] = 0;
+
+ /// three options for y:
+ /// this->starts[0] <= ny/2
+ /// ny / 2 < this->starts[0] +this->clayout->subsizes[0] < this->sizes[0] - ny/2
+ /// this->starts[0] >= this->sizes[0] - ny/2
+ /// we don't care about saving the ny/2 mode, because of symmetry
+ hsize_t y0 = this->clayout->starts[0];
+ hsize_t y1 = this->clayout->starts[0] + this->clayout->subsizes[0];
+ memshape[0] = this->clayout->subsizes[0];
+ if (y1 <= ny/2)
+ {
+ count[0] = this->clayout->subsizes[0];
+ offset[0] = y0;
+ memoffset[0] = 0;
+ }
+ else
+ {
+ if (y0 < ny/2)
+ {
+ count[0] = ny/2 - y0;
+ offset[0] = y0;
+ memoffset[0] = 0;
+ }
+ else
+ {
+ if (y1 <= this->clayout->sizes[0] - ny/2 + 1)
+ { // y0 < y1 therefore y0 <= this->clayout->sizes[0] - ny/2
+ count[0] = 0;
+ offset[0] = ny/2;
+ memoffset[0] = 0;
+ }
+ else
+ {
+ if (y0 <= this->clayout->sizes[0] - ny/2)
+ {
+ count[0] = y1 - (this->clayout->sizes[0] - ny/2);
+ offset[0] = ny/2;
+ memoffset[0] = this->clayout->subsizes[0] - count[0];
+ }
+ else
+ {
+ count[0] = this->clayout->subsizes[0];
+ offset[0] = y0;
+ memoffset[0] = 0;
+ }
+ }
+ }
+ }
+ DEBUG_MSG("count[0] = %ld, offset[0] = %ld\n",
+ count[0], offset[0]);
+ /// for z, we need to take into account that there are
+ /// both positive and negative modes
+ for (int cz = 0; cz < 2; cz++)
+ {
+ count [1] = nz/2;
+ offset[1] = cz*nz/2;
+ memshape [1] = this->clayout->sizes[1];
+ memoffset[1] = cz*(this->clayout->sizes[1] - nz/2);
+ DEBUG_MSG("cz = %d, count[1] + offset[1] = %ld\n",
+ cz, count[1] + offset[1]);
+
+ //now write data
+ mspace = H5Screate_simple(ndim(fc), memshape, NULL);
+ H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
+ H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+ H5Dwrite(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
+ H5Sclose(mspace);
+ }
+
+
+ /* close file data space */
+ H5Sclose(fspace);
+ /* close data set */
+ H5Dclose(dset_id);
+ /* close file */
+ H5Fclose(file_id);
+ return EXIT_SUCCESS;
+}
+
template <typename rnumber,
field_backend be,
diff --git a/bfps/cpp/field.hpp b/bfps/cpp/field.hpp
index a52d2a56f2527070568e54cf971b5d072520b8a3..fe6f8975e63ad1acd766b26baa84aca20ef92cc2 100644
--- a/bfps/cpp/field.hpp
+++ b/bfps/cpp/field.hpp
@@ -103,6 +103,13 @@ class field
const hid_t group,
const std::string field_name,
const int iteration);
+ int write_filtered(
+ const std::string fname,
+ const std::string field_name,
+ const int iteration,
+ const int nx,
+ const int ny,
+ const int nz);
int io_binary(
const std::string fname,