Rocky Solver SDK 2024 R2

rocky_sph_interaction

Last update: 17.07.2025
1#pragma once
2
3// Includes =======================================================================================
4#include <rocky20/api/device/api_backend.hpp>
5#include <rocky20/api/sph/rocky_sph_element.hpp>
6#include <rocky20/device/device_model.hpp>
7#include <rocky20/device/device_triangles.hpp>
8#include <rocky20/sph/device_sph.hpp>
9#include <rocky20/sph/sph_data.h>
10
11namespace device {
12__host__ __device__ inline void update_sph_element_acceleration(
13 real_sph3* element_entry, const real_sph3& acceleration)
14{
15 backend::atomic_add(&(element_entry->x), acceleration.x);
16 backend::atomic_add(&(element_entry->y), acceleration.y);
17 backend::atomic_add(&(element_entry->z), acceleration.z);
18}
19} // namespace device
20
21// ================================================================================================
22// RockySPHInteraction
23// ================================================================================================
24
25template <class NearType>
26struct RockySPHInteraction
27{
28 // Interacting relative data ------------------------------------------------------------------
29 inline ROCKY_FUNCTIONS float3 get_home_velocity() const
30 {
31 const auto near_periodic_position = [this] {
32 if constexpr (std::is_same_v<NearType, cuda_particle_sph>)
33 return this->near_element->periodic_position;
34 else
35 return PeriodicPosition {};
36 }();
37
38 float3 velocity_copy = { this->home_element->ux,
39 this->home_element->uy,
40 this->home_element->uz };
41 const auto& periodic_domain = rocky_model->gen->periodic_domain;
42 if (periodic_domain.enabled() && periodic_domain.is_cylindrical())
43 periodic_domain.get_cylindric_periodic_position(
44 velocity_copy.x,
45 velocity_copy.y,
46 velocity_copy.z,
47 this->cnt_per,
48 this->home_element->periodic_position,
49 near_periodic_position,
50 false);
51 return velocity_copy;
52 }
53
54 inline ROCKY_FUNCTIONS float3 get_distance_vector() const
55 {
56 if constexpr (std::is_same_v<NearType, cuda_particle_sph>)
57 {
58 // Distance between SPH elements
59 const auto& periodic_domain = rocky_model->gen->periodic_domain;
60 float3 home_coord { home_element->get_adjusted_coordinates(current_time) };
61 if (periodic_domain.enabled())
62 {
63 periodic_domain.get_periodic_position(
64 home_coord.x,
65 home_coord.y,
66 home_coord.z,
67 cnt_per,
68 home_element->periodic_position,
69 near_element->periodic_position,
70 false);
71 }
72 return home_coord - near_element->get_adjusted_coordinates(current_time);
73 }
74 else
75 return {};
76 }
77
78 inline ROCKY_FUNCTIONS float get_distance() const
79 {
80 return rocky::vector::get_norm(this->get_distance_vector());
81 }
82
83 // Interacting data ---------------------------------------------------------------------------
84 inline ROCKY_FUNCTIONS RockySPHElementDevice get_home_element() const
85 {
86 RockySPHElementDevice i_home_element(this->rocky_model, this->sph_model);
87 i_home_element.reset(this->home_index, this->sph_model->prt_sph + this->home_index);
88 return i_home_element;
89 }
90
91 // Contact data (resettable)
92 const cuda_particle_sph* home_element;
93 const NearType* near_element;
94 int home_index, near_index;
95 PeriodicPosition cnt_per;
96
97 // Setup data (non-resettable)
98 SDeviceModel* rocky_model;
99 SPHDeviceModel* sph_model;
100 double current_time;
101
102 ROCKY_FUNCTIONS explicit RockySPHInteraction(
103 SDeviceModel* rocky_model, SPHDeviceModel* sph_model, double current_time)
104 : rocky_model(rocky_model)
105 , sph_model(sph_model)
106 , current_time(current_time)
107 {
108 }
109};
110
111// ================================================================================================
112// RockySPHElementInteraction
113// ================================================================================================
114
115struct RockySPHElementInteraction : RockySPHInteraction<cuda_particle_sph>
116{
117 // Interacting relative data ------------------------------------------------------------------
118 // Velocity difference between the SPH elements in contact
119 inline ROCKY_FUNCTIONS float3 calculate_elements_relative_velocity()
120 {
121 const auto& periodic_domain = rocky_model->gen->periodic_domain;
122 float3 home_velocity = { home_element->ux, home_element->uy, home_element->uz };
123 if (periodic_domain.is_cylindrical())
124 periodic_domain.get_cylindric_periodic_position(
125 home_velocity.x,
126 home_velocity.y,
127 home_velocity.z,
128 cnt_per,
129 home_element->periodic_position,
130 near_element->periodic_position,
131 false);
132 return home_velocity - float3 { near_element->ux, near_element->uy, near_element->uz };
133 }
134
135 inline ROCKY_FUNCTIONS float3 get_near_velocity() const
136 {
137 float3 velocity_copy = { this->near_element->ux,
138 this->near_element->uy,
139 this->near_element->uz };
140 const auto& periodic_domain = rocky_model->gen->periodic_domain;
141 if (periodic_domain.enabled() && periodic_domain.is_cylindrical())
142 periodic_domain.get_cylindric_periodic_position(
143 velocity_copy.x,
144 velocity_copy.y,
145 velocity_copy.z,
146 this->cnt_per,
147 this->near_element->periodic_position,
148 this->home_element->periodic_position,
149 false);
150 return velocity_copy;
151 }
152
153 // Interacting data ---------------------------------------------------------------------------
154 inline ROCKY_FUNCTIONS RockySPHElementDevice get_near_element() const
155 {
156 RockySPHElementDevice i_near_element(this->rocky_model, this->sph_model);
157 i_near_element.reset(this->near_index, this->sph_model->prt_sph + this->near_index);
158 return i_near_element;
159 }
160
161 // Forces & Accelerations ---------------------------------------------------------------------
162 inline ROCKY_FUNCTIONS void add_force(float3 force)
163 {
164 // For SPH elements, the force vector actually holds acceleration in those elements.
165 // So, for SPH elements, we must extract acceleration from force.
166 //TODO Refinement
167 force /= this->sph_model->gen_sph->get_mass();
168 this->add_acceleration(force);
169 }
170
171 inline ROCKY_FUNCTIONS void add_acceleration(float3 acceleration)
172 {
173 // For SPH elements, the force vector actually holds acceleration in those elements.
174 if (this->home_element->enabled() && !this->home_element->frozen())
175 {
176 const auto& periodic_domain = rocky_model->gen->periodic_domain;
177 if (periodic_domain.is_cylindrical())
178 periodic_domain.get_cylindric_periodic_position(
179 acceleration.x,
180 acceleration.y,
181 acceleration.z,
182 this->cnt_per,
183 home_element->periodic_position,
184 near_element->periodic_position,
185 true);
186 device::update_sph_element_acceleration(sph_model->frc_sph + home_index, acceleration);
187 }
188
189 if (this->near_element->enabled() && !this->near_element->frozen())
190 device::update_sph_element_acceleration(sph_model->frc_sph + near_index, -acceleration);
191 }
192
193 inline ROCKY_FUNCTIONS void reset_home_element(int home_index)
194 {
195 // Home
196 this->home_index = home_index;
197 this->home_element = this->sph_model->prt_sph + this->home_index;
198 }
199
200 inline ROCKY_FUNCTIONS void reset_near_element(int near_index, PeriodicPosition cnt_per)
201 {
202 // Near
203 this->near_index = near_index;
204 this->near_element = this->sph_model->prt_sph + this->near_index;
205
206 // Periodic stuff
207 this->cnt_per = cnt_per;
208 }
209
210 ROCKY_FUNCTIONS explicit RockySPHElementInteraction(
211 SDeviceModel* rocky_model, SPHDeviceModel* sph_model, double current_time)
212 : RockySPHInteraction<cuda_particle_sph>(rocky_model, sph_model, current_time)
213 {
214 }
215};
216
217// ================================================================================================
218// RockySPHTriangleInteraction
219// ================================================================================================
220
221struct RockySPHTriangleInteraction : RockySPHInteraction<cuda_triangle_old>
222{
223 // Interacting relative data ------------------------------------------------------------------
224 inline ROCKY_FUNCTIONS float3 get_near_velocity() const
225 {
226 return this->boundary_velocity;
227 }
228
229 inline ROCKY_FUNCTIONS float3 get_unit_vector() const
230 {
231 return this->unit_vector;
232 }
233
234 inline ROCKY_FUNCTIONS float3 get_nearest_point_to_home() const
235 {
236 return this->inside_point;
237 }
238
239 inline ROCKY_FUNCTIONS float get_distance() const
240 {
241 return this->distance;
242 }
243
244 // Forces & accelerations ---------------------------------------------------------------------
245 inline ROCKY_FUNCTIONS void add_force(float3 force)
246 {
247 // For SPH elements, the force vector actually holds acceleration in those elements.
248 // Here we receive a force, which is straight applied to the boundaries triangles.
249 // For the home element, we must extract the acceleration.
250 //TODO Refinement
251 float3 acceleration = force / this->sph_model->gen_sph->get_mass();
252 this->add_force_and_acceleration(force, acceleration);
253 }
254
255 inline ROCKY_FUNCTIONS void add_acceleration(float3 acceleration)
256 {
257 // For SPH elements, the force vector actually holds acceleration in those elements.
258 // Here we receive an acceleration, which is straight applied to the home element.
259 // For the boundaries triangles, we must extract the force.
260 //TODO Refinement
261 float3 force = acceleration * this->sph_model->gen_sph->get_mass();
262 this->add_force_and_acceleration(force, acceleration);
263 }
264
265 inline ROCKY_FUNCTIONS void add_force_and_acceleration(
266 float3 home_acceleration, float3 near_force)
267 {
268 RockyBoundaryTriangle triangle(this->rocky_model, this->near_index);
269
270 // reduction factor for SPH elements with multiple triangle interactions
271 // (1 / number of simultaneous triangle interactions)
272 const float reduction = 1.0f / max(1.0f, this->sph_model->bnd_force_coefficient[this->home_index]);
273 home_acceleration *= reduction;
274 near_force *= reduction;
275
276 // Home element acceleration
277 if (this->home_element->enabled() && !this->home_element->frozen())
278 {
279 const auto& periodic_domain = rocky_model->gen->periodic_domain;
280 if (periodic_domain.is_cylindrical())
281 periodic_domain.get_cylindric_periodic_position(
282 home_acceleration.x,
283 home_acceleration.y,
284 home_acceleration.z,
285 this->cnt_per,
286 home_element->periodic_position,
287 {},
288 true);
289 device::update_sph_element_acceleration(
290 sph_model->frc_sph + home_index, home_acceleration);
291 }
292
293 // Near triangle: calculating moment & power
294 const auto& boundary_matrix = triangle.get_boundary_matrix();
295 double3 aux = { this->inside_point.x - boundary_matrix.dx,
296 this->inside_point.y - boundary_matrix.dy,
297 this->inside_point.z - boundary_matrix.dz };
298 double3 moment = { (aux.y * near_force.z - aux.z * near_force.y),
299 (aux.z * near_force.x - aux.x * near_force.z),
300 (aux.x * near_force.y - aux.y * near_force.x) };
301
302 // General boundary data
303 auto& boundary = triangle.get_boundary();
304
305 // Freely moving boundaries
306 if (boundary.free_f)
307 {
308 backend::atomic_add(&boundary.fx_inst, (cu_real) -near_force.x);
309 backend::atomic_add(&boundary.fy_inst, (cu_real) -near_force.y);
310 backend::atomic_add(&boundary.fz_inst, (cu_real) -near_force.z);
311 }
312 if (boundary.free_m)
313 {
314 backend::atomic_add(&boundary.mx_inst, -moment.x);
315 backend::atomic_add(&boundary.my_inst, -moment.y);
316 backend::atomic_add(&boundary.mz_inst, -moment.z);
317 }
318
319 // Record boundary statistics
320 auto* data_recorder = this->sph_model->boundary_data_recorder;
321 if (data_recorder && data_recorder->enabled)
322 {
323#ifdef __CUDA_ARCH__
324 //TODO Refinement
325 data_recorder->update_forces_data<CudaBackend>(
326 this->rocky_model->boundary_scalars,
327 this->rocky_model->triangle_scalars,
328 { -near_force.x, -near_force.y, -near_force.z },
329 -moment,
330 cast_vector<cu_real3>(this->boundary_velocity),
331 triangle.get_normal(),
332 this->sph_model->gen_sph->dt_sph,
333 triangle.get_boundary_index(),
334 this->near_index);
335#else
336 //TODO Refinement
337 data_recorder->update_forces_data<OMPBackend>(
338 this->rocky_model->boundary_scalars,
339 this->rocky_model->triangle_scalars,
340 { -near_force.x, -near_force.y, -near_force.z },
341 -moment,
342 cast_vector<cu_real3>(this->boundary_velocity),
343 triangle.get_normal(),
344 this->sph_model->gen_sph->dt_sph,
345 triangle.get_boundary_index(),
346 this->near_index);
347#endif
348 }
349 }
350
351 // Triangle Information (calculated)
352 real_sph distance;
353 real_sph normal_velocity, tangential_velocity;
354 real_sph3 inside_point, unit_vector, boundary_velocity;
355
356 ROCKY_FUNCTIONS explicit RockySPHTriangleInteraction(
357 SDeviceModel* rocky_model, SPHDeviceModel* sph_model, double current_time)
358 : RockySPHInteraction<cuda_triangle_old>(rocky_model, sph_model, current_time)
359 {
360 }
361
362 inline ROCKY_FUNCTIONS void reset_triangle_contact(
363 cuda_sph_contact* contact, PeriodicPosition cnt_per)
364 {
365 // Home
366 this->home_index = contact->nh;
367 this->home_element = this->sph_model->prt_sph + this->home_index;
368
369 // Near
370 this->near_index = contact->nn;
371 RockyBoundaryTriangle triangle(this->rocky_model, contact->nn);
372 this->near_element = triangle.triangle;
373
374 // Periodic stuff
375 this->cnt_per = cnt_per;
376
377 this->calculate_triangle_distance_information();
378 }
379
380 // Near triangle distance information
381 inline ROCKY_FUNCTIONS void calculate_triangle_distance_information()
382 {
383 const RockyBoundaryTriangle triangle(this->rocky_model, this->near_index);
384 real_sph3 position_copy = { this->home_element->x,
385 this->home_element->y,
386 this->home_element->z };
387 const auto& periodic_domain = rocky_model->gen->periodic_domain;
388 if (periodic_domain.enabled())
389 periodic_domain.get_periodic_position(
390 position_copy.x,
391 position_copy.y,
392 position_copy.z,
393 this->cnt_per,
394 home_element->periodic_position,
395 {},
396 false);
397
398 // Distance, inside point and normal unit vector
399 this->distance = device::get_triangle_orthogonal_distance(
400 triangle, position_copy, this->inside_point, this->unit_vector);
401
402 // Triangle velocity
403 cu_real3 triangle_velocity = (rocky_model->gen->flags.has_boundaries_prescribed_velocity())
404 ? rocky_model->triangle_scalars->get_known_scalar<double3>(
405 rocky20::KnownTriangleScalars::Velocity, this->near_index)
406 : double3 { 0.0, 0.0, 0.0 };
407
408 // Boundary velocity
409 this->boundary_velocity = cast_vector<real_sph3>(device::get_boundary_velocity_at_contact(
410 triangle.get_boundary(),
411 triangle.get_boundary_matrix(),
412 triangle_velocity,
413 cast_vector<double3>(this->inside_point),
414 this->current_time));
415
416 // Tangential velocity
417 real_sph3 relative_velocity = { this->boundary_velocity.x - this->home_element->ux,
418 this->boundary_velocity.y - this->home_element->uy,
419 this->boundary_velocity.z - this->home_element->uz };
420 this->normal_velocity = dot(relative_velocity, this->unit_vector);
421 this->tangential_velocity = get_norm(relative_velocity - this->normal_velocity * this->unit_vector);
422 }
423};

Connect with Ansys