Merge branch 'feature/write_slab' into 'master'

Feature/write slab

See merge request !148

cpp/field.hpp

@@ -89,6 +89,11 @@ class field
                 const std::string field_name,
                 const int iteration,
                 const bool read = true);
+        int write_rslab(
+                const std::string fname,
+                const std::string field_name,
+                const int iteration,
+                const int thickness);
         int io_database(
                 const std::string fname,
                 const std::string field_name,

cpp/field_io.cpp

@@ -345,6 +345,152 @@ int field<rnumber, be, fc>::io(
     return EXIT_SUCCESS;
 }
 
+ /** Writes a slab of a real field to a file. The slab starts at index 0 
+  * and ends at thickness in the lowest dimension (z-direction). 
+  * All other dimensions are kept untouched.
+ */
+template <typename rnumber,
+          field_backend be,
+          field_components fc>
+int field<rnumber, be, fc>::write_rslab(
+        const std::string fname,
+        const std::string field_name,
+        const int iteration,
+        const int thickness)
+{
+    /* file dataset has same dimensions as field */
+    TIMEZONE("field::write_slab");
+    assert(this->real_space_representation);
+    hid_t file_id, dset_id, plist_id;
+    file_id = H5I_BADID;
+    dset_id = H5I_BADID;
+    plist_id = H5I_BADID;
+    std::string dset_name = (
+            "/slab/" + field_name +
+            "/" + std::to_string(iteration));
+
+    /* open/create file */
+    start_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+    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);
+    if (file_id < 0) {
+        DEBUG_MSG("couldn't open file");
+        throw std::runtime_error(
+                std::string("Couldn't open file for field I/O.\n")
+                + std::string("file name = ") + fname + "\n");
+    }
+    H5Pclose(plist_id);
+    finish_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+
+    /* generic space initialization */
+    hid_t fspace, mspace;
+    hsize_t count[ndim(fc)], offset[ndim(fc)], dims[ndim(fc)];
+    hsize_t memoffset[ndim(fc)], memshape[ndim(fc)];
+
+    for (unsigned int i=0; i<ndim(fc); i++)
+    {
+        count[i] = this->rlayout->subsizes[i];
+        offset[i] = this->rlayout->starts[i];
+        dims[i] = this->rlayout->sizes[i];
+        memshape[i] = this->rmemlayout->subsizes[i];
+        memoffset[i] = 0;
+    }
+    dims[0] = thickness; //Here the thickness of the slab is set
+    if (offset[0]>=dims[0])
+    {
+        memshape[0] = 0;
+        count[0] = 0;
+    }
+    else 
+    {
+        if (offset[0]+count[0]>=dims[0])
+        {
+            memshape[0] = dims[0] - offset[0];
+            count[0] = dims[0] - offset[0];
+        }
+    }
+
+    start_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+    mspace = H5Screate_simple(ndim(fc), memshape, NULL);
+    H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
+
+    /* open/create data set */
+    if (!H5Lexists(file_id, "slab", H5P_DEFAULT))
+    {
+        hid_t gid_tmp = H5Gcreate(
+                file_id, "slab",
+                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+        H5Gclose(gid_tmp);
+    }
+
+    if (!H5Lexists(file_id, ("slab/" + field_name).c_str(), H5P_DEFAULT))
+    {
+        hid_t gid_tmp = H5Gcreate(
+                file_id, ("slab/" + field_name).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->real_space_representation ? this->rnumber_H5T : this->cnumber_H5T),
+                fspace,
+                H5P_DEFAULT,
+                H5P_DEFAULT,
+                H5P_DEFAULT);
+        assert(dset_id > 0);
+    }
+    finish_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+    /* both dset_id and fspace should now have sane values */
+
+    start_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+    /* both dset_id and fspace should now have sane values */
+    /* check file space */
+    int ndims_fspace = H5Sget_simple_extent_dims(fspace, dims, NULL);
+    variable_used_only_in_assert(ndims_fspace);
+    assert(((unsigned int)(ndims_fspace)) == ndim(fc));
+    for (unsigned int i=0; i<ndim(fc); i++)
+    {
+        offset[i] = this->rlayout->starts[i];
+        assert(dims[i] == this->rlayout->sizes[i]);
+    }
+    H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+    H5Dwrite(dset_id, this->rnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
+
+    H5Sclose(mspace);
+    H5Sclose(fspace);
+    /* close data set */
+    H5Dclose(dset_id);
+    /* ensure all processes are finished writing before exiting the method */
+    MPI_Barrier(this->comm);
+    finish_mpi_profiling_zone(turtle_mpi_pcontrol::HDF5);
+    /* close file */
+    H5Fclose(file_id);
+    return EXIT_SUCCESS;
+}
+
 template <typename rnumber,
           field_backend be,
           field_components fc>