## On this page # 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};