particle sample output arrays have wrong shape
The output datasets of sample_from_particles_system and sample_particles_system_position should have the same shape as the corresponding datasets from the checkpoint file.
At the moment they have an extra dimension of size 1.
Activity
- Author Maintainer
Suggested handling of particle dataset shape:
- C++ code reads initial condition of shape (..., d), and initial rhs of shape (n, ..., d) (where n is the order of the multi-step Adams Bashforth scheme).
- C++ code asserts that the (...) is the same for both.
- C++ code generates particle state array with shape (N, d), where N is the total number of particles (product of values in (...)). similar rhs array.
- C++ code does whatever computations.
- C++ saves checkpoints, particle datasets are shaped (..., d) and (n, ..., d).
- C++ saves samples, datasets are shaped (...) for scalar data, (..., 3) for vector data and (..., 3, 3) for tensor data.
What is needed: the shape (...) needs to be saved in the particle_system object (or something like that), and it should be accessible to
particles_outputandparticles_output_sampling. Ideally, the shape should be set when the output objects are initialized, since the shape will not change once a DNS is started. - Maintainer
I am not sure to get it. currently (...) is the number of particles right? so what you like to have instead? a vector?
- Author Maintainer
(...) could be something like (m1, m2, ..., mk), and N = m1m2...mk, with arbitrary k.
So shapes:
- state (m1, m2, ..., mk, d)
- rhs (n, m1, m2, ..., mk, d)
- scalar data (m1, m2, ..., mk) (field template parameter
ONE) - vector data (m1, m2, ..., mk, 3) (field template parameter
THREE) - tensor data (m1, m2, ..., mk, 3, 3) (field template parameter
THREExTHREE).
In particular, please look at
NSVEparticlesP2P::do_stats, and how I usesave_datasetthere.Edited by Cristian Lalescu - Maintainer
Ok, I see. You will use this information only for IO? so you will get it when you load and use it when you write, that's all?
- Author Maintainer
yes, only for IO.
- Maintainer
I am currently looking to it. But, it is a little difficult. Either it can end in unbalanced number of particles per process, or some process will have to access the data multiple times.
For example, let
(m1, m2, ..., mk)be(10,3,5,8)(so N = 1200). If P = 11, we should split around 108 particles per process. How can we load 108 particles? it will need multiple loads right? Or maybe there is something that I did not see? - Author Maintainer
I am looking at
particles_input_hdf5.hpp, lines 149 through 173. I will post an update of this in 10 mins. - Maintainer
I have these modifications:
{ TIMEZONE("state"); hid_t dset = H5Dopen(particle_file, inDatanameState.c_str(), H5P_DEFAULT); assert(dset >= 0); hid_t dspace = H5Dget_space(dset); // copy? assert(dspace >= 0); hid_t space_dim = H5Sget_simple_extent_ndims(dspace); assert(space_dim >= 2); std::vector<hsize_t> state_dim_array(space_dim); int hdfret = H5Sget_simple_extent_dims(dspace, &state_dim_array[0], NULL); assert(hdfret >= 0); // Last value is the position dim of the particles assert(state_dim_array.back() == size_particle_positions); // The data shape contains the dims except the last one dataShape = state_dim_array; dataShape.pop_back(); nb_total_particles = 1; for (size_t idx_dim = 0; idx_dim < state_dim_array.size()-1; ++idx_dim){ nb_total_particles *= state_dim_array[idx_dim]; } hdfret = H5Sclose(dspace); assert(hdfret >= 0); hdfret = H5Dclose(dset); assert(hdfret >= 0); } { TIMEZONE("rhs"); hid_t dset = H5Dopen(particle_file, inDatanameRhs.c_str(), H5P_DEFAULT); assert(dset >= 0); hid_t dspace = H5Dget_space(dset); // copy? assert(dspace >= 0); hid_t rhs_dim = H5Sget_simple_extent_ndims(dspace); assert(rhs_dim >= 0); std::vector<hsize_t> rhs_dim_array(rhs_dim); int hdfret = H5Sget_simple_extent_dims(dspace, &rhs_dim_array[0], NULL); assert(hdfret >= 0); assert(rhs_dim_array.back() == size_particle_rhs); assert(rhs_dim_array.front() >= 1); assert(rhs_dim == 2 + dataShape.size()); for(size_t idx_shape = 0 ; idx_shape < dataShape.size() ; ++idx_shape){ assert(dataShape[idx_shape] == rhs_dim_array[idx_shape+1]); } nb_rhs = rhs_dim_array[0]; hdfret = H5Sclose(dspace); assert(hdfret >= 0); hdfret = H5Dclose(dset); assert(hdfret >= 0); } particles_utils::IntervalSplitter<hsize_t> load_splitter(nb_total_particles, nb_processes, my_rank); static_assert(std::is_same<real_number, double>::value || std::is_same<real_number, float>::value, "real_number must be double or float"); const hid_t type_id = (sizeof(real_number) == 8?H5T_NATIVE_DOUBLE:H5T_NATIVE_FLOAT); /// Load the data std::unique_ptr<real_number[]> split_particles_positions(new real_number[load_splitter.getMySize()*size_particle_positions]); { TIMEZONE("state-read"); hid_t dset = H5Dopen(particle_file, inDatanameState.c_str(), H5P_DEFAULT); assert(dset >= 0); hid_t rspace = H5Dget_space(dset); assert(rspace >= 0); std::vector<hsize_t> offset = get_offset_with_shape(load_splitter.getMyOffset(), dataShape); offset.push_back(0); std::vector<hsize_t> mem_dims = get_size_with_shape(load_splitter.getMySize(), dataShape, offset); mem_dims.push_back(size_particle_positions); hid_t mspace = H5Screate_simple(mem_dims.size(), &mem_dims[0], NULL); assert(mspace >= 0); int rethdf = H5Sselect_hyperslab(rspace, H5S_SELECT_SET, offset, NULL, mem_dims, NULL); assert(rethdf >= 0); rethdf = H5Dread(dset, type_id, mspace, rspace, H5P_DEFAULT, split_particles_positions.get()); assert(rethdf >= 0); rethdf = H5Sclose(rspace); assert(rethdf >= 0); rethdf = H5Dclose(dset); assert(rethdf >= 0); } std::vector<std::unique_ptr<real_number[]>> split_particles_rhs(nb_rhs); { TIMEZONE("rhs-read"); hid_t dset = H5Dopen(particle_file, inDatanameRhs.c_str(), H5P_DEFAULT); assert(dset >= 0); std::vector<hsize_t> offset = get_offset_with_shape(load_splitter.getMyOffset(), dataShape); offset.push_front(0); offset.push_back(0); std::vector<hsize_t> mem_dims = get_size_with_shape(load_splitter.getMySize(), dataShape, offset); mem_dims.push_front(1); mem_dims.push_back(size_particle_positions); for(hsize_t idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){ hid_t rspace = H5Dget_space(dset); assert(rspace >= 0); split_particles_rhs[idx_rhs].reset(new real_number[load_splitter.getMySize()*size_particle_rhs]); offset[0] = idx_rhs; hid_t mspace = H5Screate_simple( mem_dims.size(), &mem_dims[0], NULL); assert(mspace >= 0); int rethdf = H5Sselect_hyperslab( rspace, H5S_SELECT_SET, offset, NULL, mem_dims, NULL); assert(rethdf >= 0); rethdf = H5Dread(dset, type_id, mspace, rspace, H5P_DEFAULT, split_particles_rhs[idx_rhs].get()); assert(rethdf >= 0); rethdf = H5Sclose(mspace); assert(rethdf >= 0); rethdf = H5Sclose(rspace); assert(rethdf >= 0); } int rethdf = H5Dclose(dset); assert(rethdf >= 0); } - Maintainer
And I was on:
std::vector<hsize_t> offset = get_offset_with_shape(load_splitter.getMyOffset(), dataShape); offset.push_back(0); std::vector<hsize_t> mem_dims = get_size_with_shape(load_splitter.getMySize(), dataShape, offset); mem_dims.push_back(size_particle_positions);And I think it will not be possible in one read or it will be unbalanced.
- Author Maintainer
Here's my replacement of lines 149 through 173:
{ TIMEZONE("state-read"); hid_t dset = H5Dopen(particle_file, inDatanameState.c_str(), H5P_DEFAULT); assert(dset >= 0); hid_t rspace = H5Dget_space(dset); assert(rspace >= 0); // get number of dimensions int rspace_ndims = H5Sget_simple_extent_ndims(rspace); // get sizes hsize_t *rspace_dims = new hsize_t[rspace_ndims]; H5Sget_simple_extent_dims(rspace, rspace_dims, NULL); hsize_t *rspace_offset = new hsize_t[rspace_ndims]; // for state component, offset is always 0 rspace_offset[0] = 0; hsize_t total_nparticle_offset = load_splitter.getMyOffset(); for (int dim_counter = 1; dim_counter < rspace_ndims; dim_counter++) { rspace_offset[dim_counter] = total_particle_offset % rspace_dims[dim_counter]; total_particle_offset = total_particle_offset / rspace_dims[dim_counter]; } hsize_t mem_dims[2] = {load_splitter.getMySize(), size_particle_positions}; hid_t mspace = H5Screate_simple(2, &mem_dims[0], NULL); assert(mspace >= 0); int rethdf = H5Sselect_hyperslab(rspace, H5S_SELECT_SET, rspace_offset, NULL, mem_dims, NULL); assert(rethdf >= 0); rethdf = H5Dread(dset, type_id, mspace, rspace, H5P_DEFAULT, split_particles_positions.get()); assert(rethdf >= 0); rethdf = H5Sclose(rspace); assert(rethdf >= 0); rethdf = H5Dclose(dset); assert(rethdf >= 0); delete[] rspace_dims; delete[] rspace_offset; }Edited by Cristian Lalescu - Maintainer
But this is not over here, right?
hsize_t mem_dims[2] = {load_splitter.getMySize(), size_particle_positions}; hid_t mspace = H5Screate_simple(2, &mem_dims[0], NULL); - Author Maintainer
mspace and mem_dims stay the same. HDF5 will not complain, because the sizes will be compatible (unless I did something wrong with the computation of the offset).
- Author Maintainer
Hm. My code is incomplete, you're right. In select_hyperslab something other than mem_dims should be used. I'll post the correction in a few minutes.
- Author Maintainer
Here's the corrected code:
{ TIMEZONE("state-read"); hid_t dset = H5Dopen(particle_file, inDatanameState.c_str(), H5P_DEFAULT); assert(dset >= 0); hid_t rspace = H5Dget_space(dset); assert(rspace >= 0); // get number of dimensions int rspace_ndims = H5Sget_simple_extent_ndims(rspace); // get sizes hsize_t *rspace_dims = new hsize_t[rspace_ndims]; H5Sget_simple_extent_dims(rspace, rspace_dims, NULL); hsize_t *rspace_offset = new hsize_t[rspace_ndims]; hsize_t *rspace_count = new hsize_t[rspace_ndims]; // for state component, offset is always 0 rspace_offset[0] = 0; rspace_count[0] = size_particle_positions; hsize_t total_nparticle_offset = load_splitter.getMyOffset(); hsize_t total_nparticle_count = load_splitter.getMySize(); for (int dim_counter = 1; dim_counter < rspace_ndims; dim_counter++) { rspace_offset[dim_counter] = total_particle_offset % rspace_dims[dim_counter]; total_particle_offset = total_particle_offset / rspace_dims[dim_counter]; rspace_count[dim_counter] = total_particle_count % rspace_dims[dim_counter]; total_particle_count = total_particle_count / rspace_dims[dim_counter]; } hsize_t mem_dims[2] = {load_splitter.getMySize(), size_particle_positions}; hid_t mspace = H5Screate_simple(2, &mem_dims[0], NULL); assert(mspace >= 0); int rethdf = H5Sselect_hyperslab(rspace, H5S_SELECT_SET, rspace_offset, NULL, rspace_count, NULL); assert(rethdf >= 0); rethdf = H5Dread(dset, type_id, mspace, rspace, H5P_DEFAULT, split_particles_positions.get()); assert(rethdf >= 0); rethdf = H5Sclose(rspace); assert(rethdf >= 0); rethdf = H5Dclose(dset); assert(rethdf >= 0); delete[] rspace_dims; delete[] rspace_offset; }Edited by Cristian Lalescu - Author Maintainer
As far as I can tell, my code selects contiguous chunks of the dataset, and then it puts them into differently shaped array of the same size. HDF5 should not complain: see 7.4.1.4. of http://davis.lbl.gov/Manuals/HDF5-1.8.7/UG/12_Dataspaces.html
Edited by Cristian Lalescu - Author Maintainer
I just noticed. Why do you use
&mem_dims[0]instead ofmem_dims? - Maintainer
Ok, I see that was what I was missing, that is great, so let me try to update my version. I use
&mem_dims[0]because it is a vector not a pointer (I could use .data()). - Maintainer
Ha but in the previous code there is not reason.
- Maintainer
Do we agree that the loop:
for (int dim_counter = 1; dim_counter < rspace_ndims; dim_counter++)should be in the other order? from
rspace_ndims-1to1? - Author Maintainer
I'm not sure.
If I have 135, and I want to write it in decimal:
I compute 135 % 10, this gives me 5 I divide 135 by 10, this gives me 13 2) I compute 13 % 10, this gives me 3 I divide 13 by 10, this gives me 1 3) I compute 1 % 10, this gives me 1
Aren't these the digits we want to compute?