diff --git a/src/core/visualization/paraview/vtk_diffusion_grid.cc b/src/core/visualization/paraview/vtk_diffusion_grid.cc index 8e816877f..dd7c918fa 100644 --- a/src/core/visualization/paraview/vtk_diffusion_grid.cc +++ b/src/core/visualization/paraview/vtk_diffusion_grid.cc @@ -42,8 +42,10 @@ VtkDiffusionGrid::VtkDiffusionGrid(const std::string& name, if (param->export_visualization) { auto* tinfo = ThreadInfo::GetInstance(); data_.resize(tinfo->GetMaxThreads()); + piece_boxes_z_.resize(tinfo->GetMaxThreads()); } else { data_.resize(1); + piece_boxes_z_.resize(1); } #pragma omp parallel for schedule(static, 1) @@ -104,6 +106,10 @@ VtkDiffusionGrid::~VtkDiffusionGrid() { bool VtkDiffusionGrid::IsUsed() const { return used_; } // ----------------------------------------------------------------------------- +// ToDo (tobias): +// * remove debug code (print statements) +// * remove now unnecessary branching +// * add comments void VtkDiffusionGrid::Update(const DiffusionGrid* grid) { used_ = true; @@ -147,41 +153,60 @@ void VtkDiffusionGrid::Update(const DiffusionGrid* grid) { } #pragma omp parallel for schedule(static, 1) - for (uint64_t i = 0; i < num_pieces_; ++i) { + for (uint64_t i = 0; i < piece_boxes_z_.size(); ++i) { uint64_t piece_elements; auto* e = piece_extents_[i].data(); - if (i < num_pieces_ - 1) { - piece_elements = piece_boxes_z_ * xy_num_boxes; - data_[i]->SetDimensions(num_boxes[0], num_boxes[1], piece_boxes_z_); + if (i < piece_boxes_z_.size() - 1) { + piece_elements = piece_boxes_z_[i] * xy_num_boxes; + data_[i]->SetDimensions(num_boxes[0], num_boxes[1], piece_boxes_z_[i]); data_[i]->SetExtent(e[0], e[1], e[2], e[3], e[4], - e[4] + piece_boxes_z_ - 1); + e[4] + piece_boxes_z_[i] - 1); } else { - piece_elements = piece_boxes_z_last_ * xy_num_boxes; - data_[i]->SetDimensions(num_boxes[0], num_boxes[1], piece_boxes_z_last_); + piece_elements = piece_boxes_z_[i] * xy_num_boxes; + data_[i]->SetDimensions(num_boxes[0], num_boxes[1], piece_boxes_z_[i]); data_[i]->SetExtent(e[0], e[1], e[2], e[3], e[4], - e[4] + piece_boxes_z_last_ - 1); + e[4] + piece_boxes_z_[i] - 1); } - real_t piece_origin_z = origin_z + box_length * piece_boxes_z_ * i; + // Compute partial sum of boxes until index i + auto sum = + std::accumulate(piece_boxes_z_.begin(), piece_boxes_z_.begin() + i, 0); + real_t piece_origin_z = origin_z + box_length * sum; data_[i]->SetOrigin(origin_x, origin_y, piece_origin_z); data_[i]->SetSpacing(box_length, box_length, box_length); + // // Debug + // std::cout << "sum: " << sum << std::endl; + // std::cout << "piece_origin_z : " << piece_origin_z << std::endl; if (concentration_array_idx_ != -1) { auto* co_ptr = const_cast(grid->GetAllConcentrations()); auto elements = static_cast(piece_elements); auto* array = static_cast( data_[i]->GetPointData()->GetArray(concentration_array_idx_)); - if (i < num_pieces_ - 1) { - array->SetArray(co_ptr + (elements * i), elements, 1); + auto offset = sum * xy_num_boxes; + if (i < piece_boxes_z_.size() - 1) { + // // Debug + // std::cout << "Offset : " << offset << std::endl; + // std::cout << "Value (begin) : " << co_ptr[offset] << std::endl; + // std::cout << "Value (end): " << co_ptr[offset + elements - 1] + // << std::endl; + array->SetArray(co_ptr + offset, elements, 1); } else { - array->SetArray(co_ptr + total_boxes - elements, elements, 1); + // // Debug + // std::cout << "Offset : " << offset << std::endl; + // std::cout << "Value (begin) : " << co_ptr[offset] << std::endl; + // std::cout << "Value (end): " << co_ptr[offset + elements - 1] + // << std::endl; + array->SetArray(co_ptr + offset, elements, 1); } + // // Debug + // std::cout << "elements : " << elements << std::endl; } if (gradient_array_idx_ != -1) { auto gr_ptr = const_cast(grid->GetAllGradients()); auto elements = static_cast(piece_elements * 3); auto* array = static_cast( data_[i]->GetPointData()->GetArray(gradient_array_idx_)); - if (i < num_pieces_ - 1) { + if (i < piece_boxes_z_.size() - 1) { array->SetArray(gr_ptr + (elements * i), elements, 1); } else { array->SetArray(gr_ptr + total_boxes - elements, elements, 1); @@ -196,47 +221,66 @@ void VtkDiffusionGrid::WriteToFile(uint64_t step) const { auto filename_prefix = Concat(name_, "-", step); ParallelVtiWriter writer; - writer(sim->GetOutputDir(), filename_prefix, data_, num_pieces_, + writer(sim->GetOutputDir(), filename_prefix, data_, piece_boxes_z_.size(), whole_extent_, piece_extents_); } // ----------------------------------------------------------------------------- void VtkDiffusionGrid::Dissect(uint64_t boxes_z, uint64_t num_pieces_target) { - if (num_pieces_target == 1) { - piece_boxes_z_last_ = boxes_z; - piece_boxes_z_ = 1; - num_pieces_ = 1; - } else if (boxes_z <= num_pieces_target) { - piece_boxes_z_last_ = 2; - piece_boxes_z_ = 1; - num_pieces_ = std::max(1ULL, boxes_z - 1); - } else { - auto boxes_per_piece = static_cast(boxes_z) / num_pieces_target; - piece_boxes_z_ = static_cast(std::ceil(boxes_per_piece)); - num_pieces_ = boxes_z / piece_boxes_z_; - piece_boxes_z_last_ = boxes_z - (num_pieces_ - 1) * piece_boxes_z_; + auto boxes_per_piece = static_cast(boxes_z) / num_pieces_target; + auto min_slices = static_cast(std::floor(boxes_per_piece)); + auto leftover = boxes_z % num_pieces_target; + // // Debug print info of all variables + // std::cout << "boxes_z : " << boxes_z << std::endl; + // std::cout << "num_pieces_target : " << num_pieces_target << std::endl; + // std::cout << "boxes_per_piece : " << boxes_per_piece << std::endl; + // std::cout << "min_slices : " << min_slices << std::endl; + // std::cout << "leftover : " << leftover << std::endl; + // Distribute leftover slices + for (uint64_t i = 0; i < piece_boxes_z_.size(); ++i) { + piece_boxes_z_[i] = min_slices; + if (leftover > 0) { + piece_boxes_z_[i]++; + leftover--; + } + } + // Shorten vector by removing all tailing zeros + auto it = std::find(piece_boxes_z_.begin(), piece_boxes_z_.end(), 0); + if (it != piece_boxes_z_.end()) { + piece_boxes_z_.resize(std::distance(piece_boxes_z_.begin(), it)); } - assert(num_pieces_ > 0); - assert(piece_boxes_z_last_ >= 1); - assert((num_pieces_ - 1) * piece_boxes_z_ + piece_boxes_z_last_ == boxes_z); + + // Check dissection + auto sum = std::accumulate(piece_boxes_z_.begin(), piece_boxes_z_.end(), 0); + assert(sum == boxes_z); + // // Debug: Calculate sum of all entries + // std::cout << "Sum : " << sum << std::endl; + // std::cout << "boxes_z : " << boxes_z << std::endl; + + // // Debug print the vector + // std::cout << "piece_boxes_z_ : "; + // for (auto i : piece_boxes_z_) { + // std::cout << i << " "; + // } + // std::cout << std::endl; } // ----------------------------------------------------------------------------- void VtkDiffusionGrid::CalcPieceExtents( const std::array& num_boxes) { - piece_extents_.resize(num_pieces_); - if (num_pieces_ == 1) { + piece_extents_.resize(piece_boxes_z_.size()); + if (piece_boxes_z_.size() == 1) { piece_extents_[0] = whole_extent_; return; } - int c = piece_boxes_z_; + int c = piece_boxes_z_[0]; piece_extents_[0] = {{0, static_cast(num_boxes[0]) - 1, 0, static_cast(num_boxes[1]) - 1, 0, c}}; - for (uint64_t i = 1; i < num_pieces_ - 1; ++i) { + for (uint64_t i = 1; i < piece_boxes_z_.size() - 1; ++i) { piece_extents_[i] = {{0, static_cast(num_boxes[0]) - 1, 0, static_cast(num_boxes[1]) - 1, c, - c + static_cast(piece_boxes_z_)}}; - c += piece_boxes_z_; + c + static_cast(piece_boxes_z_[i])}}; + c += piece_boxes_z_[i]; } piece_extents_.back() = {{0, static_cast(num_boxes[0]) - 1, 0, static_cast(num_boxes[1]) - 1, c, diff --git a/src/core/visualization/paraview/vtk_diffusion_grid.h b/src/core/visualization/paraview/vtk_diffusion_grid.h index d63db9558..aea223c55 100644 --- a/src/core/visualization/paraview/vtk_diffusion_grid.h +++ b/src/core/visualization/paraview/vtk_diffusion_grid.h @@ -55,9 +55,7 @@ class VtkDiffusionGrid { // vtkImageData objects for parallel processing. // ParaView calls the partition pieces. - uint64_t num_pieces_; - uint64_t piece_boxes_z_; - uint64_t piece_boxes_z_last_; + std::vector piece_boxes_z_; std::array whole_extent_; std::vector> piece_extents_;