BioDynaMo  v1.03.58-27764645
mechanical_forces_op_cuda.cc
Go to the documentation of this file.
1 // -----------------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 CERN & University of Surrey for the benefit of the
4 // BioDynaMo collaboration. All Rights Reserved.
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 //
9 // See the LICENSE file distributed with this work for details.
10 // See the NOTICE file distributed with this work for additional information
11 // regarding copyright ownership.
12 //
13 // -----------------------------------------------------------------------------
14 
15 #include <vector>
16 
18 #include "core/agent/cell.h"
24 #include "core/resource_manager.h"
25 #include "core/shape.h"
26 #include "core/simulation.h"
27 #include "core/util/log.h"
28 #include "core/util/thread_info.h"
29 #include "core/util/timing.h"
30 #include "core/util/type.h"
31 #include "core/util/vtune_helper.h"
32 
33 namespace bdm {
34 
35 // -----------------------------------------------------------------------------
36 void IsNonSphericalObjectPresent(const Agent* agent, bool* answer) {
37  if (agent->GetShape() != Shape::kSphere) {
38  *answer = true;
39  }
40 }
41 
42 namespace detail {
43 
44 // -----------------------------------------------------------------------------
45 struct InitializeGPUData : public Functor<void, Agent*, AgentHandle> {
47 
48  real_t* cell_movements = nullptr;
49  real_t* cell_positions = nullptr;
50  real_t* cell_diameters = nullptr;
51  real_t* cell_adherence = nullptr;
53  uint32_t* cell_boxid = nullptr;
54  real_t* mass = nullptr;
55  uint32_t* successors = nullptr;
56 
57  std::vector<AgentHandle::ElementIdx_t> offset;
58 
59  uint32_t* starts = nullptr;
60  uint16_t* lengths = nullptr;
61  uint64_t* timestamps = nullptr;
62  uint64_t* current_timestamp = nullptr;
63  uint32_t* num_boxes_axis = nullptr;
65 
66  uint64_t allocated_num_agents = 0;
67  uint64_t allocated_num_boxes = 0;
68 
70 
71  virtual ~InitializeGPUData();
72 
73  void Initialize(uint64_t num_agents, uint64_t num_boxes,
74  const std::vector<AgentHandle::ElementIdx_t>& offs,
76 
77  void operator()(Agent* agent, AgentHandle ah) override;
78 
79  private:
80  void FreeAgentBuffers();
81 
82  void FreeGridBuffers();
83 };
84 
85 // -----------------------------------------------------------------------------
87 
88 // -----------------------------------------------------------------------------
90  if (current_timestamp != nullptr) {
93  }
94 
95  if (allocated_num_agents != 0) {
97  }
98 
99  if (allocated_num_boxes != 0) {
100  FreeGridBuffers();
101  }
102 }
103 
104 // -----------------------------------------------------------------------------
106  uint64_t num_agents, uint64_t num_boxes,
107  const std::vector<AgentHandle::ElementIdx_t>& offs,
109  if (current_timestamp == nullptr) {
112  }
113 
114  if (allocated_num_agents < num_agents) {
115  if (allocated_num_agents != 0) {
117  }
118  allocated_num_agents = num_agents * 1.25;
127  }
128 
129  if (allocated_num_boxes < num_boxes) {
130  if (allocated_num_boxes != 0) {
131  FreeGridBuffers();
132  }
133  allocated_num_boxes = num_boxes * 1.25;
137  }
138 
139  offset = offs;
140  grid = g;
141 }
142 
143 // -----------------------------------------------------------------------------
153 }
154 
155 // -----------------------------------------------------------------------------
160 }
161 
162 // -----------------------------------------------------------------------------
164  // Check if there are any non-spherical objects in our simulation, because
165  // GPU accelerations currently supports only sphere-sphere interactions
168  Log::Fatal("MechanicalForcesOpCuda",
169  "\nWe detected a non-spherical object during the GPU "
170  "execution. This is currently not supported.");
171  return;
172  }
173  auto* cell = bdm_static_cast<Cell*>(agent);
174  auto idx = offset[ah.GetNumaNode()] + ah.GetElementIdx();
175  auto idxt3 = idx * 3;
176  mass[idx] = cell->GetMass();
177  cell_diameters[idx] = cell->GetDiameter();
178  cell_adherence[idx] = cell->GetAdherence();
179  const auto& tf = cell->GetTractorForce();
180  const auto& pos = cell->GetPosition();
181 
182  for (uint64_t i = 0; i < 3; ++i) {
183  cell_tractor_force[idxt3 + i] = tf[i];
184  cell_positions[idxt3 + i] = pos[i];
185  }
186 
187  cell_boxid[idx] = cell->GetBoxIdx();
188 
189  // populate successor list
190  auto nid = ah.GetNumaNode();
191  auto el = ah.GetElementIdx();
192  successors[idx] = offset[grid->successors_.data_[nid][el].GetNumaNode()] +
193  grid->successors_.data_[nid][el].GetElementIdx();
194 }
195 
196 } // namespace detail
197 
198 // -----------------------------------------------------------------------------
199 struct UpdateCPUResults : public Functor<void, Agent*, AgentHandle> {
200  real_t* cell_movements = nullptr;
201  std::vector<AgentHandle::ElementIdx_t> offset;
202 
204  const std::vector<AgentHandle::ElementIdx_t>& offs);
205  virtual ~UpdateCPUResults();
206 
207  void operator()(Agent* agent, AgentHandle ah) override;
208 };
209 
210 // -----------------------------------------------------------------------------
212  real_t* cm, const std::vector<AgentHandle::ElementIdx_t>& offs) {
213  cell_movements = cm;
214  offset = offs;
215 }
216 
217 // -----------------------------------------------------------------------------
219 
220 // -----------------------------------------------------------------------------
222  auto* param = Simulation::GetActive()->GetParam();
223  auto* cell = bdm_static_cast<Cell*>(agent);
224  auto idx = offset[ah.GetNumaNode()] + ah.GetElementIdx();
225  idx *= 3;
226  Real3 new_pos = {cell_movements[idx], cell_movements[idx + 1],
227  cell_movements[idx + 2]};
228  cell->UpdatePosition(new_pos);
229  if (param->bound_space) {
230  ApplyBoundingBox(agent, param->bound_space, param->min_bound,
231  param->max_bound);
232  }
233 }
234 
235 // -----------------------------------------------------------------------------
237  // Timing timer("MechanicalForcesOpCuda::SetUp");
238  auto* sim = Simulation::GetActive();
239  auto* grid = dynamic_cast<UniformGridEnvironment*>(sim->GetEnvironment());
240  auto* rm = sim->GetResourceManager();
241 
242  if (!grid) {
243  Log::Fatal(
244  "MechanicalForcesOpCuda::operator()",
245  "MechanicalForcesOpCuda only works with UniformGridEnvironement.");
246  }
247 
248  auto num_numa_nodes = ThreadInfo::GetInstance()->GetNumaNodes();
249  std::vector<AgentHandle::ElementIdx_t> offset(num_numa_nodes);
250  offset[0] = 0;
251  for (int nn = 1; nn < num_numa_nodes; nn++) {
252  offset[nn] = offset[nn - 1] + rm->GetNumAgents(nn - 1);
253  }
254 
255  auto total_num_agents = rm->GetNumAgents();
256  auto num_boxes = grid->boxes_.size();
257 
258  if (i_ == nullptr) {
260  }
261  i_->Initialize(total_num_agents, num_boxes, offset, grid);
262  {
263  // Timing timer("MechanicalForcesOpCuda::toColumnar");
264  rm->ForEachAgentParallel(1000, *i_);
265  }
266 
267  *i_->current_timestamp = grid->timestamp_;
268 #pragma omp parallel for
269  for (uint64_t i = 0; i < grid->boxes_.size(); ++i) {
270  auto& box = grid->boxes_[i];
271  i_->timestamps[i] = box.timestamp_;
272  if (box.timestamp_ == *i_->current_timestamp) {
273  i_->lengths[i] = box.length_;
274  i_->starts[i] =
275  offset[box.start_.GetNumaNode()] + box.start_.GetElementIdx();
276  }
277  }
278  grid->GetNumBoxesAxis(i_->num_boxes_axis);
279 }
280 
281 // -----------------------------------------------------------------------------
283  auto* sim = Simulation::GetActive();
284  auto* grid = dynamic_cast<UniformGridEnvironment*>(sim->GetEnvironment());
285  auto* param = sim->GetParam();
286  auto* rm = sim->GetResourceManager();
287 
288  uint32_t total_num_agents = rm->GetNumAgents();
289  auto num_boxes = grid->boxes_.size();
290  // If this is the first time we perform physics on GPU using CUDA
291  if (cdo_ == nullptr) {
292  // Allocate 25% more memory so we don't need to reallocate GPU memory
293  // for every (small) change
294  total_num_agents_ = static_cast<uint32_t>(1.25 * total_num_agents);
295  num_boxes_ = static_cast<uint32_t>(1.25 * num_boxes);
296 
297  // Allocate required GPU memory
299  } else {
300  // If the number of agents increased
301  if (total_num_agents >= total_num_agents_) {
302  Log::Info("MechanicalForcesOpCuda",
303  "\nThe number of cells increased signficantly (from ",
304  total_num_agents_, " to ", total_num_agents,
305  "), agent we allocate bigger GPU buffers\n");
306  total_num_agents_ = static_cast<uint32_t>(1.25 * total_num_agents);
308  }
309 
310  // If the neighbor grid size increased
311  if (num_boxes >= num_boxes_) {
312  Log::Info("MechanicalForcesOpCuda",
313  "\nThe number of boxes increased signficantly (from ",
314  num_boxes_, " to ", "), so we allocate bigger GPU buffers\n");
315  num_boxes_ = static_cast<uint32_t>(1.25 * num_boxes);
317  }
318  }
319 
320  real_t squared_radius =
321  grid->GetLargestAgentSize() * grid->GetLargestAgentSize();
322 
323  // Timing timer("MechanicalForcesOpCuda::Kernel");
326  i_->cell_adherence, i_->cell_boxid, i_->mass, param->simulation_time_step,
327  param->simulation_max_displacement, squared_radius, total_num_agents,
330 }
331 
332 // -----------------------------------------------------------------------------
334  cdo_->Sync();
335  // Timing timer("MechanicalForcesOpCuda::TearDown");
338 }
339 
340 } // namespace bdm
timing.h
shape.h
bdm::MechanicalForcesOpCudaKernel
Definition: mechanical_forces_op_cuda_kernel.h:29
bdm::detail::InitializeGPUData::~InitializeGPUData
virtual ~InitializeGPUData()
Definition: mechanical_forces_op_cuda.cc:89
bdm::UpdateCPUResults::operator()
void operator()(Agent *agent, AgentHandle ah) override
Definition: mechanical_forces_op_cuda.cc:221
bdm::ThreadInfo::GetInstance
static ThreadInfo * GetInstance()
Definition: thread_info.cc:21
bdm::detail::InitializeGPUData::Initialize
void Initialize(uint64_t num_agents, uint64_t num_boxes, const std::vector< AgentHandle::ElementIdx_t > &offs, UniformGridEnvironment *g)
Definition: mechanical_forces_op_cuda.cc:105
agent_handle.h
bdm::MechanicalForcesOpCudaKernel::ResizeCellBuffers
void ResizeCellBuffers(uint32_t num_cells)
Definition: mechanical_forces_op_cuda_kernel.h:82
bdm
Definition: agent.cc:39
bdm::detail::InitializeGPUData::FreeGridBuffers
void FreeGridBuffers()
Definition: mechanical_forces_op_cuda.cc:156
bdm::MechanicalForcesOpCudaKernel::ResizeGridBuffers
void ResizeGridBuffers(uint32_t num_boxes)
Definition: mechanical_forces_op_cuda_kernel.h:84
bdm::detail::InitializeGPUData::current_timestamp
uint64_t * current_timestamp
Definition: mechanical_forces_op_cuda.cc:62
thread_info.h
bdm::UpdateCPUResults::offset
std::vector< AgentHandle::ElementIdx_t > offset
Definition: mechanical_forces_op_cuda.cc:201
bdm::ThreadInfo::GetNumaNodes
int GetNumaNodes() const
Returns the number of NUMA nodes on this machine.
Definition: thread_info.h:48
vtune_helper.h
bdm::MechanicalForcesOpCuda::SetUp
void SetUp() override
Definition: mechanical_forces_op_cuda.cc:236
bdm::UniformGridEnvironment
A class that represents Cartesian 3D grid.
Definition: uniform_grid_environment.h:58
bdm::real_t
double real_t
Definition: real_t.h:21
bdm::CudaAllocPinned
void CudaAllocPinned(T **d, uint64_t elements)
Definition: cuda_pinned_memory.h:33
bdm::UniformGridEnvironment::boxes_
ParallelResizeVector< Box > boxes_
Definition: uniform_grid_environment.h:651
bdm::UpdateCPUResults
Definition: mechanical_forces_op_cuda.cc:199
bdm::detail::InitializeGPUData::offset
std::vector< AgentHandle::ElementIdx_t > offset
Definition: mechanical_forces_op_cuda.cc:57
mechanical_forces_op_cuda.h
type.h
bdm::Agent::GetShape
virtual Shape GetShape() const =0
bdm::UniformGridEnvironment::successors_
AgentVector< AgentHandle > successors_
Definition: uniform_grid_environment.h:676
bdm::Agent
Contains code required by all agents.
Definition: agent.h:79
bdm::ApplyBoundingBox
void ApplyBoundingBox(Agent *agent, Param::BoundSpaceMode mode, real_t lb, real_t rb)
Definition: bound_space_op.h:26
uniform_grid_environment.h
bdm::UpdateCPUResults::~UpdateCPUResults
virtual ~UpdateCPUResults()
bdm::MechanicalForcesOpCuda::cdo_
MechanicalForcesOpCudaKernel * cdo_
Definition: mechanical_forces_op_cuda.h:42
bdm::Functor
Definition: functor.h:24
bdm::detail::InitializeGPUData::successors
uint32_t * successors
Definition: mechanical_forces_op_cuda.cc:55
bdm::detail::InitializeGPUData::cell_movements
real_t * cell_movements
Definition: mechanical_forces_op_cuda.cc:48
bdm::detail::InitializeGPUData::timestamps
uint64_t * timestamps
Definition: mechanical_forces_op_cuda.cc:61
bdm::Log::Info
static void Info(const std::string &location, const Args &... parts)
Prints information message.
Definition: log.h:55
bdm::Simulation::GetResourceManager
ResourceManager * GetResourceManager()
Returns the ResourceManager instance.
Definition: simulation.cc:261
bdm::detail::InitializeGPUData::mass
real_t * mass
Definition: mechanical_forces_op_cuda.cc:54
bdm::MechanicalForcesOpCudaKernel::Sync
void Sync() const
Definition: mechanical_forces_op_cuda_kernel.h:81
bdm::kSphere
@ kSphere
Definition: shape.h:20
bdm::CudaFreePinned
void CudaFreePinned(void *p)
Definition: cuda_pinned_memory.h:35
bdm::detail::InitializeGPUData::allocated_num_agents
uint64_t allocated_num_agents
Definition: mechanical_forces_op_cuda.cc:66
bdm::detail::InitializeGPUData::is_non_spherical_object
bool is_non_spherical_object
Definition: mechanical_forces_op_cuda.cc:46
bdm::Log::Fatal
static void Fatal(const std::string &location, const Args &... parts)
Prints fatal error message.
Definition: log.h:115
bdm::UpdateCPUResults::cell_movements
real_t * cell_movements
Definition: mechanical_forces_op_cuda.cc:200
bound_space_op.h
bdm::MechanicalForcesOpCuda::total_num_agents_
uint32_t total_num_agents_
Definition: mechanical_forces_op_cuda.h:45
bdm::detail::InitializeGPUData::FreeAgentBuffers
void FreeAgentBuffers()
Definition: mechanical_forces_op_cuda.cc:144
log.h
bdm::detail::InitializeGPUData::cell_positions
real_t * cell_positions
Definition: mechanical_forces_op_cuda.cc:49
bdm::detail::InitializeGPUData::num_boxes_axis
uint32_t * num_boxes_axis
Definition: mechanical_forces_op_cuda.cc:63
bdm::detail::InitializeGPUData::cell_tractor_force
real_t * cell_tractor_force
Definition: mechanical_forces_op_cuda.cc:52
bdm::MechanicalForcesOpCuda::operator()
void operator()() override
Definition: mechanical_forces_op_cuda.cc:282
bdm::MechanicalForcesOpCuda::num_boxes_
uint32_t num_boxes_
Definition: mechanical_forces_op_cuda.h:44
bdm::detail::InitializeGPUData::InitializeGPUData
InitializeGPUData()
bdm::detail::InitializeGPUData::allocated_num_boxes
uint64_t allocated_num_boxes
Definition: mechanical_forces_op_cuda.cc:67
bdm::detail::InitializeGPUData::grid
UniformGridEnvironment * grid
Definition: mechanical_forces_op_cuda.cc:64
environment.h
bdm::ResourceManager::ForEachAgentParallel
virtual void ForEachAgentParallel(Functor< void, Agent * > &function, Functor< bool, Agent * > *filter=nullptr)
Definition: resource_manager.cc:92
simulation.h
bdm::Simulation::GetParam
const Param * GetParam() const
Returns the simulation parameters.
Definition: simulation.cc:271
bdm::MechanicalForcesOpCudaKernel::LaunchMechanicalForcesKernel
void LaunchMechanicalForcesKernel(const real_t *positions, const real_t *diameter, const real_t *tractor_force, const real_t *adherence, const uint32_t *box_id, const real_t *mass, const real_t timestep, const real_t max_displacement, const real_t squared_radius, const uint32_t num_agents, uint32_t *starts, uint16_t *lengths, uint64_t *timestamps, uint64_t current_timestamp, uint32_t *successors, uint32_t *num_boxes_axis, real_t *cell_movements)
Definition: mechanical_forces_op_cuda_kernel.h:72
bdm::detail::InitializeGPUData::cell_adherence
real_t * cell_adherence
Definition: mechanical_forces_op_cuda.cc:51
bdm::MathArray< real_t, 3 >
bdm::detail::InitializeGPUData::cell_boxid
uint32_t * cell_boxid
Definition: mechanical_forces_op_cuda.cc:53
bdm::AgentHandle::GetNumaNode
NumaNode_t GetNumaNode() const
Definition: agent_handle.h:44
resource_manager.h
bdm::Simulation::GetActive
static Simulation * GetActive()
This function returns the currently active Simulation simulation.
Definition: simulation.cc:67
bdm::UpdateCPUResults::UpdateCPUResults
UpdateCPUResults(real_t *cm, const std::vector< AgentHandle::ElementIdx_t > &offs)
Definition: mechanical_forces_op_cuda.cc:211
bdm::IsNonSphericalObjectPresent
void IsNonSphericalObjectPresent(const Agent *agent, bool *answer)
Definition: mechanical_forces_op_cuda.cc:36
cuda_pinned_memory.h
bdm::MechanicalForcesOpCuda::i_
detail::InitializeGPUData * i_
Definition: mechanical_forces_op_cuda.h:43
bdm::detail::InitializeGPUData::operator()
void operator()(Agent *agent, AgentHandle ah) override
Definition: mechanical_forces_op_cuda.cc:163
bdm::AgentHandle::GetElementIdx
ElementIdx_t GetElementIdx() const
Definition: agent_handle.h:45
cell.h
bdm::detail::InitializeGPUData::lengths
uint16_t * lengths
Definition: mechanical_forces_op_cuda.cc:60
bdm::detail::InitializeGPUData::starts
uint32_t * starts
Definition: mechanical_forces_op_cuda.cc:59
bdm::detail::InitializeGPUData
Definition: mechanical_forces_op_cuda.cc:45
bdm::MechanicalForcesOpCuda::TearDown
void TearDown() override
Definition: mechanical_forces_op_cuda.cc:333
bdm::detail::InitializeGPUData::cell_diameters
real_t * cell_diameters
Definition: mechanical_forces_op_cuda.cc:50
bdm::AgentHandle
Definition: agent_handle.h:29