bugfix: write_filtered correct for vectors and tensors

TurTLE/test/test_write_filtered.py

@@ -24,7 +24,19 @@ def main():
     df = h5py.File('data_for_write_filtered.h5', 'r')
     bla0 = df['scal_field_full/complex/0'][:, 0]
     bla1 = df['scal_field_z0slice/complex/0'][...]
-    assert(np.max(np.abs(bla1 - bla0)) < 1e-5)
+    max_dist = np.max(np.abs(bla1 - bla0))
+    print('maximum distance for scalar is ', max_dist)
+    assert(max_dist < 1e-5)
+    bla0 = df['vec_field_full/complex/0'][:, 0]
+    bla1 = df['vec_field_z0slice/complex/0'][...]
+    max_dist = np.max(np.abs(bla1 - bla0))
+    print('maximum distance for vector is ', max_dist)
+    assert(max_dist < 1e-5)
+    bla0 = df['tens_field_full/complex/0'][:, 0]
+    bla1 = df['tens_field_z0slice/complex/0'][...]
+    max_dist = np.max(np.abs(bla1 - bla0))
+    print('maximum distance for tensor is ', max_dist)
+    assert(max_dist < 1e-5)
     print('SUCCESS!!! z=0 slice agrees between different datasets')
     return None
 

cpp/field.cpp

@@ -672,12 +672,9 @@ int field<rnumber, be, fc>::write_filtered(
     }
     // these are dimensions of dataset, needed
     // to create dataset
-    //dims[0] = nz;
     dims[0] = ny;
     dims[1] = nz;
     dims[2] = nx/2+1;
-    if (nz == 1)
-        dims[1] = nx/2+1;
 
     /* open/create data set */
     if (!H5Lexists(file_id, field_name.c_str(), H5P_DEFAULT))
@@ -702,10 +699,16 @@ int field<rnumber, be, fc>::write_filtered(
     else
     {
         if (nz == 1)
+        {
+            hsize_t temp_dims[ndim(fc)-1];
+            temp_dims[0] = dims[0];
+            for (unsigned int i=1; i<ndim(fc)-1; i++)
+                temp_dims[i] = dims[i+1];
             fspace = H5Screate_simple(
                     ndim(fc)-1,
-                    dims,
+                    temp_dims,
                     NULL);
+        }
         else
             fspace = H5Screate_simple(
                     ndim(fc),
@@ -728,7 +731,8 @@ int field<rnumber, be, fc>::write_filtered(
     {
         assert(((unsigned int)(ndims_fspace)) == ndim(fc)-1);
         assert(dims[0] == fdims[0]);
-        assert(dims[1] == fdims[1]);
+        for (unsigned int i=1; i<ndim(fc)-1; i++)
+            assert(dims[i+1] == fdims[i]);
     }
     else
     {
@@ -830,7 +834,8 @@ int field<rnumber, be, fc>::write_filtered(
         //now write data
         mspace = H5Screate_simple(ndim(fc), memshape, NULL);
         H5Sselect_hyperslab(mspace, H5S_SELECT_SET, memoffset, NULL, count, NULL);
-        count[1] = nx/2+1;
+        for (unsigned int i=1; i<ndim(fc)-1; i++)
+            count[i] = count[i+1];
         H5Sselect_hyperslab(fspace, H5S_SELECT_SET, offset, NULL, count, NULL);
         H5Dwrite(dset_id, this->cnumber_H5T, mspace, fspace, H5P_DEFAULT, this->data);
         H5Sclose(mspace);

cpp/full_code/write_filtered_test.cpp

@@ -110,6 +110,98 @@ int write_filtered_test<rnumber>::do_work(void)
 
     // deallocate
     delete scal_field;
+
+    // allocate
+    field<rnumber, FFTW, THREE> *vec_field = new field<rnumber, FFTW, THREE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            FFTW_ESTIMATE);
+
+    // fill up vec_field
+    vec_field->real_space_representation = true;
+    vec_field->RLOOP(
+            [&](ptrdiff_t rindex,
+                ptrdiff_t xindex,
+                ptrdiff_t yindex,
+                ptrdiff_t zindex){
+            vec_field->rval(rindex, 0) = rdist(rgen);
+            vec_field->rval(rindex, 1) = rdist(rgen);
+            vec_field->rval(rindex, 2) = rdist(rgen);
+            });
+    vec_field->dft();
+
+    vec_field->write_filtered(
+            "data_for_write_filtered.h5",
+            "vec_field_z0slice",
+            this->iteration,
+            this->nx,
+            this->ny,
+            1);
+
+    vec_field->write_filtered(
+            "data_for_write_filtered.h5",
+            "vec_field_full",
+            this->iteration,
+            this->nx,
+            this->ny,
+            this->nz);
+
+    vec_field->write_filtered(
+            "data_for_write_filtered.h5",
+            "vec_field_half",
+            this->iteration,
+            this->nx/2,
+            this->ny/2,
+            this->nz/2);
+
+    // deallocate
+    delete vec_field;
+
+    // allocate
+    field<rnumber, FFTW, THREExTHREE> *tens_field = new field<rnumber, FFTW, THREExTHREE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            FFTW_ESTIMATE);
+
+    // fill up tens_field
+    tens_field->real_space_representation = true;
+    tens_field->RLOOP(
+            [&](ptrdiff_t rindex,
+                ptrdiff_t xindex,
+                ptrdiff_t yindex,
+                ptrdiff_t zindex){
+            for (int cc = 0; cc<3; cc++)
+            for (int ccc = 0; ccc<3; ccc++)
+                tens_field->rval(rindex, cc, ccc) = rdist(rgen);
+            });
+    tens_field->dft();
+
+    tens_field->write_filtered(
+            "data_for_write_filtered.h5",
+            "tens_field_z0slice",
+            this->iteration,
+            this->nx,
+            this->ny,
+            1);
+
+    tens_field->write_filtered(
+            "data_for_write_filtered.h5",
+            "tens_field_full",
+            this->iteration,
+            this->nx,
+            this->ny,
+            this->nz);
+
+    tens_field->write_filtered(
+            "data_for_write_filtered.h5",
+            "tens_field_half",
+            this->iteration,
+            this->nx/2,
+            this->ny/2,
+            this->nz/2);
+
+    // deallocate
+    delete tens_field;
     return EXIT_SUCCESS;
 }