Merge branch 'feature/slices' into develop

bfps/NavierStokes.py

@@ -263,20 +263,6 @@ class NavierStokes(_fluid_particle_base):
             field_H5T = 'H5T_NATIVE_FLOAT'
         elif self.dtype == np.float64:
             field_H5T = 'H5T_NATIVE_DOUBLE'
-        self.stat_src += self.create_stat_output(
-                '/statistics/xlines/velocity',
-                'fs->rvelocity',
-                data_type = field_H5T,
-                size_setup = """
-                    count[0] = 1;
-                    count[1] = nx;
-                    count[2] = 3;
-                    """,
-                close_spaces = False)
-        self.stat_src += self.create_stat_output(
-                '/statistics/xlines/vorticity',
-                'fs->rvorticity',
-                data_type = field_H5T)
         if self.QR_stats_on:
             self.stat_src += self.create_stat_output(
                     '/statistics/moments/trS2_Q_R',
@@ -754,12 +740,14 @@ class NavierStokes(_fluid_particle_base):
             vec_stat_datasets = ['velocity', 'vorticity']
             scal_stat_datasets = []
             for k in vec_stat_datasets:
-                time_chunk = 2**20//(8*3*self.parameters['nx']) # FIXME: use proper size of self.dtype
+                time_chunk = 2**20 // (
+                        self.dtype.itemsize*3*
+                        self.parameters['nx']*self.parameters['ny'])
                 time_chunk = max(time_chunk, 1)
-                ofile.create_dataset('statistics/xlines/' + k,
-                                     (1, self.parameters['nx'], 3),
-                                     chunks = (time_chunk, self.parameters['nx'], 3),
-                                     maxshape = (None, self.parameters['nx'], 3),
+                ofile.create_dataset('statistics/0slices/' + k + '/real',
+                                     (1, self.parameters['ny'], self.parameters['nx'], 3),
+                                     chunks = (time_chunk, self.parameters['ny'], self.parameters['nx'], 3),
+                                     maxshape = (None, self.parameters['ny'], self.parameters['nx'], 3),
                                      dtype = self.dtype)
             if self.Lag_acc_stats_on:
                 vec_stat_datasets += ['Lagrangian_acceleration']

bfps/cpp/field.cpp

@@ -542,6 +542,73 @@ int field<rnumber, be, fc>::io_database(
 }
 
 
+template <typename rnumber,
+          field_backend be,
+          field_components fc>
+int field<rnumber, be, fc>::write_0slice(
+        const hid_t group,
+        const std::string field_name,
+        const int iteration)
+{
+    TIMEZONE("field::write_0slice");
+    assert(this->real_space_representation);
+    assert(fc == THREE);
+    if (this->myrank == 0)
+    {
+        hid_t dset, wspace, mspace;
+        int ndims;
+        hsize_t count[4], offset[4], dims[4];
+        offset[0] = iteration;
+        offset[1] = 0;
+        offset[2] = 0;
+        offset[3] = 0;
+        dset = H5Dopen(
+                group,
+                ("0slices/" + field_name + "/real").c_str(),
+                H5P_DEFAULT);
+        wspace = H5Dget_space(dset);
+        ndims = H5Sget_simple_extent_dims(wspace, dims, NULL);
+        // array in memory has 2 extra x points, because FFTW
+        count[0] = 1;
+        count[1] = this->rmemlayout->sizes[1];
+        count[2] = this->rmemlayout->sizes[2];
+        count[3] = 3;
+        mspace = H5Screate_simple(ndims, count, NULL);
+        // array in file should not have the extra 2 points
+        count[1] = this->rlayout->sizes[1];
+        count[2] = this->rlayout->sizes[2];
+        // select right slice in file
+        H5Sselect_hyperslab(
+            wspace,
+            H5S_SELECT_SET,
+            offset,
+            NULL,
+            count,
+            NULL);
+        offset[0] = 0;
+        // select proper regions of memory
+        H5Sselect_hyperslab(
+            mspace,
+            H5S_SELECT_SET,
+            offset,
+            NULL,
+            count,
+            NULL);
+        H5Dwrite(
+            dset,
+            this->rnumber_H5T,
+            mspace,
+            wspace,
+            H5P_DEFAULT,
+            this->data);
+        H5Dclose(dset);
+        H5Sclose(mspace);
+        H5Sclose(wspace);
+    }
+    return EXIT_SUCCESS;
+}
+
+
 template <typename rnumber,
           field_backend be,
           field_components fc>
@@ -762,6 +829,11 @@ void field<rnumber, be, fc>::compute_rspace_stats(
         H5Sclose(wspace);
         H5Sclose(mspace);
         H5Dclose(dset);
+        if (H5Lexists(group, "0slices", H5P_DEFAULT))
+            this->write_0slice(
+                    group,
+                    dset_name,
+                    toffset);
     }
     delete[] moments;
     delete[] hist;
@@ -1039,4 +1111,4 @@ template void compute_gradient<float, FFTW, ONE, THREE, SMOOTH>(
 template void compute_gradient<double, FFTW, ONE, THREE, SMOOTH>(
         kspace<FFTW, SMOOTH> *,
         field<double, FFTW, ONE> *,
-        field<double, FFTW, THREE> *);
\ No newline at end of file
+        field<double, FFTW, THREE> *);

bfps/cpp/field.hpp

@@ -85,6 +85,11 @@ class field
                 const int toffset,
                 const bool read = true);
 
+        int write_0slice(
+                const hid_t group,
+                const std::string field_name,
+                const int iteration);
+
         /* essential FFT stuff */
         void dft();
         void ift();