From e0aafef1e3fe5d5812926b8a3d51498447f2564f Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Tue, 19 Nov 2019 16:57:52 +0100
Subject: [PATCH] bugfix: write_filtered correct for vectors and tensors

---
 TurTLE/test/test_write_filtered.py    | 14 +++-
 cpp/field.cpp                         | 17 +++--
 cpp/full_code/write_filtered_test.cpp | 92 +++++++++++++++++++++++++++
 3 files changed, 116 insertions(+), 7 deletions(-)

diff --git a/TurTLE/test/test_write_filtered.py b/TurTLE/test/test_write_filtered.py
index 6a484907..b43fda32 100644
--- a/TurTLE/test/test_write_filtered.py
+++ b/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
 
diff --git a/cpp/field.cpp b/cpp/field.cpp
index 7303d332..04de2b0c 100644
--- a/cpp/field.cpp
+++ b/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);
diff --git a/cpp/full_code/write_filtered_test.cpp b/cpp/full_code/write_filtered_test.cpp
index fa6e05c3..3e270794 100644
--- a/cpp/full_code/write_filtered_test.cpp
+++ b/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;
 }
 
-- 
GitLab