fix writing of particle state

bfps/cpp/distributed_particles.cpp

@@ -387,116 +387,46 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
             this->comm);
     delete[] yy;
     if (this->myrank == 0)
-    {
-        std::string temp_string = (std::string(this->name) +
-                                   std::string("/") +
-                                   std::string(dset_name));
-        hid_t dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT);
-        hid_t mspace, wspace;
-        hsize_t count[3], offset[3];
-        wspace = H5Dget_space(dset);
-        H5Sget_simple_extent_dims(wspace, count, NULL);
-        assert(count[2] == 3);
-        count[0] = 1;
-        offset[0] = this->iteration / this->traj_skip;
-        offset[1] = 0;
-        offset[2] = 0;
-        mspace = H5Screate_simple(3, count, NULL);
-        H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
-        H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, data);
-        H5Sclose(mspace);
-        H5Sclose(wspace);
-        H5Dclose(dset);
-    }
+        this->write_point3D_chunk(dset_name, 0, data);
     delete[] data;
 }
 
 template <int particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
         const hid_t data_file_id,
-        const char *dset_name,
-        std::unordered_map<int, single_particle_state<particle_type>> &y)
+        const bool write_rhs)
+{
+    double *temp0 = new double[this->chunk_size*this->ncomponents];
+    double *temp1 = new double[this->chunk_size*this->ncomponents];
+    for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
     {
-    //double *data = new double[this->nparticles*this->ncomponents];
-    //double *yy = new double[this->nparticles*this->ncomponents];
-    //std::fill_n(yy, this->nparticles*this->ncomponents, 0);
-    //for (int p=0; p<this->nparticles; p++)
-    //{
-    //    auto pp = y.find(p);
-    //    if (pp != y.end())
-    //        std::copy(pp->second.data,
-    //                  pp->second.data + this->ncomponents,
-    //                  yy + pp->first*this->ncomponents);
-    //}
-    //MPI_Allreduce(
-    //        yy,
-    //        data,
-    //        this->ncomponents*this->nparticles,
-    //        MPI_DOUBLE,
-    //        MPI_SUM,
-    //        this->comm);
-    //delete[] yy;
-    //if (this->myrank == 0)
-    //{
-    //    std::string temp_string = (std::string(this->name) +
-    //                               std::string("/") +
-    //                               std::string(dset_name));
-    //    hid_t dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT);
-    //    hid_t mspace, wspace;
-    //    hsize_t count[3], offset[3];
-    //    wspace = H5Dget_space(dset);
-    //    H5Sget_simple_extent_dims(wspace, count, NULL);
-    //    count[0] = 1;
-    //    offset[0] = this->iteration / this->traj_skip;
-    //    offset[1] = 0;
-    //    offset[2] = 0;
-    //    mspace = H5Screate_simple(3, count, NULL);
-    //    H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
-    //    H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, data);
-    //    H5Sclose(mspace);
-    //    H5Sclose(wspace);
-    //    H5Dclose(dset);
-    //}
-    //delete[] data;
+        //write state
+        std::fill_n(temp0, this->ncomponents*this->chunk_size, 0);
+        for (int p=0; p<this->chunk_size; p++)
+        {
+            auto pp = this->state.find(p + cindex*this->chunk_size);
+            if (pp != this->state.end())
+                std::copy(pp->second.data,
+                          pp->second.data + this->ncomponents,
+                          temp0 + pp->first*this->ncomponents);
         }
-
-template <int particle_type, class rnumber, int interp_neighbours>
-void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
-        const hid_t data_file_id,
-        const bool write_rhs)
+        MPI_Allreduce(
+                temp0,
+                temp1,
+                this->ncomponents*this->chunk_size,
+                MPI_DOUBLE,
+                MPI_SUM,
+                this->comm);
+        if (this->myrank == 0)
+            this->write_state_chunk(cindex, temp1);
+        //write rhs
+        if (this->iteration > 0)
+            for (int i=0; i<this->integration_steps; i++)
             {
-    this->write(data_file_id, "state", this->state);
-    //if (this->myrank == 0)
-    //{
-    //    if (write_rhs)
-    //    {
-    //        std::string temp_string = (
-    //                std::string("/") +
-    //                std::string(this->name) +
-    //                std::string("/rhs"));
-    //        hid_t dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT);
-    //        hid_t wspace = H5Dget_space(dset);
-    //        hsize_t count[4], offset[4];
-    //        H5Sget_simple_extent_dims(wspace, count, NULL);
-    //        //writing to last available position
-    //        offset[0] = count[0] - 1;
-    //        offset[1] = 0;
-    //        offset[2] = 0;
-    //        offset[3] = 0;
-    //        count[0] = 1;
-    //        count[1] = 1;
-    //        hid_t mspace = H5Screate_simple(4, count, NULL);
-    //        for (int i=0; i<this->integration_steps; i++)
-    //        {
-    //            offset[1] = i;
-    //            H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
-    //            H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, this->rhs[i]);
-    //        }
-    //        H5Sclose(mspace);
-    //        H5Sclose(wspace);
-    //        H5Dclose(dset);
-    //    }
-    //}
+            }
+    }
+    delete[] temp0;
+    delete[] temp1;
 }