1 // ============================================================================= 2 // PROJECT CHRONO - http://projectchrono.org 3 // 4 // Copyright (c) 2019 projectchrono.org 5 // All rights reserved. 6 // 7 // Use of this source code is governed by a BSD-style license that can be found 8 // in the LICENSE file at the top level of the distribution and at 9 // http://projectchrono.org/license-chrono.txt. 10 // 11 // ============================================================================= 12 // Authors: Conlain Kelly, Nic Olsen, Dan Negrut, Radu Serban 13 // ============================================================================= 14 15 #pragma once 16 17 #include <cstddef> 18 #include <functional> 19 #include <string> 20 #include <vector> 21 #include <algorithm> 22 #include <cmath> 23 24 #include "chrono_gpu/ChApiGpu.h" 25 #include "chrono_gpu/ChGpuDefines.h" 26 #include "chrono_gpu/physics/ChGpuBoundaryConditions.h" 27 #include "chrono_gpu/cuda/ChGpuCUDAalloc.hpp" 28 29 typedef unsigned char not_stupid_bool; 30 31 namespace chrono { 32 namespace gpu { 33 34 /// <summary> 35 /// ChSolverStateData contains information that pertains the solver, at a certain point in time. It contains 36 /// information about the current sim time, current time step, max number of spheres in each SD, etc. 37 /// This is the type of information that changes (or not) from time step to time step. 38 /// </summary> 39 class ChSolverStateData { 40 private: 41 /// variable stores at each time step the largest number of spheres in any domain; this is useful when deciding 42 /// on the number of threads in a block for a kernel call. It is also useful to settle once and for all whether 43 /// we exceed the threshold MAX_COUNT_OF_SPHERES_PER_SD. Avoids having these checks in the kernel. Stored in 44 /// managed memory. 45 unsigned int* pMaxNumberSpheresInAnySD; 46 unsigned int crntMaxNumberSpheresInAnySD; // redundant with info above, but info above is on the device 47 unsigned int largestMaxNumberSpheresInAnySD_thusFar; 48 49 /// vector of unsigned int that lives on the device; used by CUB or by anybody else that needs scrap space. 50 /// Please pay attention to the type the vector stores. 51 std::vector<char, cudallocator<char>> deviceScratchSpace; 52 53 /// current integration time step 54 float crntStepSize_SU; // DN: needs to be brought here from GranParams 55 float crntSimTime_SU; // DN: needs to be brought here from GranParams 56 public: ChSolverStateData()57 ChSolverStateData() { 58 cudaMallocManaged(&pMaxNumberSpheresInAnySD, sizeof(unsigned int)); 59 largestMaxNumberSpheresInAnySD_thusFar = 0; 60 } ~ChSolverStateData()61 ~ChSolverStateData() { cudaFree(pMaxNumberSpheresInAnySD); } pMM_maxNumberSpheresInAnySD()62 inline unsigned int* pMM_maxNumberSpheresInAnySD() { 63 return pMaxNumberSpheresInAnySD; ///< returns pointer to managed memory 64 } 65 /// keep track of the largest number of spheres that touched an SD thus far into the simulation recordCrntMaxNumberSpheresInAnySD(unsigned int someVal)66 inline void recordCrntMaxNumberSpheresInAnySD(unsigned int someVal) { 67 crntMaxNumberSpheresInAnySD = someVal; 68 if (someVal > largestMaxNumberSpheresInAnySD_thusFar) 69 largestMaxNumberSpheresInAnySD_thusFar = someVal; 70 } 71 /// reports the largest number of spheres that touched any SD thus far into the simulation MaxNumberSpheresInAnySDThusFar()72 inline unsigned int MaxNumberSpheresInAnySDThusFar() const { return largestMaxNumberSpheresInAnySD_thusFar; } 73 /// reports the largest number of spheres that touched any SD at the most recent broadphase CD function call currentMaxNumberSpheresInAnySD()74 inline unsigned int currentMaxNumberSpheresInAnySD() const { return crntMaxNumberSpheresInAnySD; } 75 76 /// return raw pointer to swath of device memory that is at least "sizeNeeded" large pDeviceMemoryScratchSpace(size_t sizeNeeded)77 inline char* pDeviceMemoryScratchSpace(size_t sizeNeeded) { 78 if (deviceScratchSpace.size() < sizeNeeded) { 79 deviceScratchSpace.resize(sizeNeeded, 0); 80 } 81 return deviceScratchSpace.data(); 82 } 83 }; 84 85 // Underlying implementation of the Chrono::Gpu system. 86 // used to control and dispatch the GPU sphere-only solver. 87 class ChSystemGpu_impl { 88 public: 89 virtual ~ChSystemGpu_impl(); 90 91 protected: 92 /// Structure with simulation parameters for sphere-based granular dynamics. 93 /// This structure is stored in CUDA unified memory so that it can be accessed from both host and device. 94 struct GranParams { 95 float stepSize_SU; ///< Timestep in SU 96 97 CHGPU_FRICTION_MODE friction_mode; ///< Which friction mode is active for the simulation 98 CHGPU_ROLLING_MODE rolling_mode; ///< Which rolling resistance model is active 99 CHGPU_TIME_INTEGRATOR time_integrator; ///< Which time integrator is active 100 101 /// Ratio of normal force to peak tangent force, also arctan(theta) where theta is the friction angle 102 /// sphere-to-sphere 103 float static_friction_coeff_s2s; 104 /// Ratio of normal force to peak tangent force, also arctan(theta) where theta is the friction angle 105 /// sphere-to-wall 106 float static_friction_coeff_s2w; 107 108 float rolling_coeff_s2s_SU; ///< Coefficient of rolling resistance sphere-to-sphere 109 float rolling_coeff_s2w_SU; ///< Coefficient of rolling resistance sphere-to-wall 110 111 float spinning_coeff_s2s_SU; ///< Coefficient of spinning resistance sphere-to-sphere 112 float spinning_coeff_s2w_SU; ///< Coefficient of spinning resistance sphere-to-wall 113 114 float Gamma_n_s2s_SU; ///< sphere-to-sphere normal contact damping coefficient, expressed in SU 115 float Gamma_n_s2w_SU; ///< sphere-to-wall normal contact damping coefficient, expressed in SU 116 float Gamma_t_s2s_SU; ///< sphere-to-sphere tangent contact damping coefficient, expressed in SU 117 float Gamma_t_s2w_SU; ///< sphere-to-wall tangent contact damping coefficient, expressed in SU 118 119 float K_n_s2s_SU; ///< sphere-to-sphere normal contact stiffness, expressed in SU 120 float K_n_s2w_SU; ///< sphere-to-wall normal contact stiffness, expressed in SU 121 float K_t_s2s_SU; ///< sphere-to-sphere tangent contact stiffness, expressed in SU 122 float K_t_s2w_SU; ///< sphere-to-wall tangent contact stiffness, expressed in SU 123 124 /// use material-based property 125 bool use_mat_based = false; 126 127 /// effective sphere-to-sphere youngs modulus, expressed in SU 128 float E_eff_s2s_SU; 129 /// effective sphere-to-wall youngs modulus, expressed in SU 130 float E_eff_s2w_SU; 131 132 /// effective sphere-to-sphere shear modulus, expressed in SU 133 float G_eff_s2s_SU; 134 /// effective sphere-to-wall shear modulus, expressed in SU 135 float G_eff_s2w_SU; 136 137 /// effective sphere-to-sphere coefficient of restitution, expressed in SU 138 float COR_s2s_SU; 139 /// effective sphere-to-wall coefficient of restitution, expressed in SU 140 float COR_s2w_SU; 141 142 unsigned int sphereRadius_SU; ///< Radius of the sphere, expressed in SU 143 float sphereInertia_by_r; ///< Moment of inertia of a sphere, normalized by the radius, expressed in SU 144 145 unsigned int SD_size_X_SU; ///< X-dimension of each subdomain box, expressed in SU 146 unsigned int SD_size_Y_SU; ///< Y-dimension of each subdomain box, expressed in SU 147 unsigned int SD_size_Z_SU; ///< Z-dimension of each subdomain box, expressed in SU 148 149 unsigned int nSpheres; ///< Total number of spheres in system, used for boundary condition multistep friction 150 unsigned int nSDs; ///< Total number of subdomains 151 152 unsigned int nSDs_X; ///< X-dimension of the big domain box in multiples of subdomains 153 unsigned int nSDs_Y; ///< Y-dimension of the big domain box in multiples of subdomains 154 unsigned int nSDs_Z; ///< Z-dimension of the big domain box in multiples of subdomains 155 156 int64_t max_x_pos_unsigned; ///< Maximum X dimension in the big domain frame ((int64_t)SD_size_X_SU * nSDs_X) 157 int64_t max_y_pos_unsigned; ///< Maximum Y dimension in the big domain frame ((int64_t)SD_size_Y_SU * nSDs_Y) 158 int64_t max_z_pos_unsigned; ///< Maximum Z dimension in the big domain frame ((int64_t)SD_size_Z_SU * nSDs_Z) 159 160 float gravAcc_X_SU; ///< X gravity in SU 161 float gravAcc_Y_SU; ///< Y gravity in SU 162 float gravAcc_Z_SU; ///< Z gravity in SU 163 164 int64_t BD_frame_X; ///< The bottom-left corner xPos of the big domain 165 int64_t BD_frame_Y; ///< The bottom-left corner yPos of the big domain 166 int64_t BD_frame_Z; ///< The bottom-left corner zPos of the big domain 167 168 int64_t BD_offset_X; /// X offset of big domain from its original frame; allows the subdomain to move 169 int64_t BD_offset_Y; /// Y offset of big domain from its original frame; allows the subdomain to move 170 int64_t BD_offset_Z; /// Z offset of big domain from its original frame; allows the subdomain to move 171 172 float cohesionAcc_s2s; ///< Constant acceleration of sphere-to-sphere cohesion 173 float adhesionAcc_s2w; ///< Accleration of adhesion 174 175 double LENGTH_UNIT; ///< 1 / C_L. Any length expressed in SU is a multiple of LENGTH_UNIT 176 double TIME_UNIT; ///< 1 / C_T. Any time quantity in SU is measured as a positive multiple of TIME_UNIT 177 double MASS_UNIT; ///< 1 / C_M. Any mass quantity is measured as a positive multiple of MASS_UNIT 178 179 /// Make clear that the underlying assumption is unit SU mass 180 constexpr static float sphere_mass_SU = 1.f; 181 182 /// Used as a safety check to determine whether a system has lost stability 183 float max_safe_vel = (float)UINT_MAX; 184 185 bool recording_contactInfo = false; ///< recording contact info 186 }; 187 188 /// Structure of pointers to kinematic quantities of the ChSystemGpu_impl. 189 /// These pointers must be in device-accessible memory. 190 struct SphereData { 191 int* sphere_local_pos_X; ///< X position relative to owner subdomain in unified memory 192 int* sphere_local_pos_Y; ///< Y position relative to owner subdomain in unified memory 193 int* sphere_local_pos_Z; ///< Z position relative to owner subdomain in unified memory 194 195 float* pos_X_dt; ///< X velocity in unified memory 196 float* pos_Y_dt; ///< Y velocity in unified memory 197 float* pos_Z_dt; ///< Z velocity in unified memory 198 199 float* sphere_Omega_X; ///< X angular velocity in unified memory. Only used if friction is present 200 float* sphere_Omega_Y; ///< Y angular velocity in unified memory. Only used if friction is present 201 float* sphere_Omega_Z; ///< Z angular velocity in unified memory. Only used if friction is present 202 203 float* sphere_acc_X; ///< X angular acceleration in unified memory 204 float* sphere_acc_Y; ///< Y angular acceleration in unified memory 205 float* sphere_acc_Z; ///< Z angular acceleration in unified memory 206 207 float* sphere_ang_acc_X; ///< X angular acceleration in unified memory 208 float* sphere_ang_acc_Y; ///< Y angular acceleration in unified memory 209 float* sphere_ang_acc_Z; ///< Z angular acceleration in unified memory 210 211 float* sphere_acc_X_old; ///< Previous step X acceleration for multistep integrators 212 float* sphere_acc_Y_old; ///< Previous step Y acceleration for multistep integrators 213 float* sphere_acc_Z_old; ///< Previous step Z acceleration for multistep integrators 214 215 float* sphere_ang_acc_X_old; ///< Previous step X angular acceleration for multistep integrators 216 float* sphere_ang_acc_Y_old; ///< Previous step Y angular acceleration for multistep integrators 217 float* sphere_ang_acc_Z_old; ///< Previous step Z angular acceleration for multistep integrators 218 219 not_stupid_bool* sphere_fixed; ///< Flags indicating whether or not a sphere is fixed 220 221 float* sphere_stats_buffer; ///< A buffer array that can store any quantity that the user wish to reduce 222 223 unsigned int* contact_partners_map; ///< Contact partners for each sphere. Only in frictional simulations 224 not_stupid_bool* contact_active_map; ///< Whether the frictional contact at an index is active 225 float3* contact_history_map; ///< Tangential history for a given contact pair. Only for multistep friction 226 float* contact_duration; ///< Duration of persistent contact between pairs 227 228 float3* normal_contact_force; ///< Track normal contact force 229 float3* tangential_friction_force; ///< Track sliding friction force 230 float3* rolling_friction_torque; ///< Track rolling friction force 231 float* char_collision_time; ///< Track characteristic collision time 232 float3* v_rot_array; ///< Track v_rot array 233 234 unsigned int* SD_NumSpheresTouching; ///< Number of particles touching each subdomain 235 unsigned int* SD_SphereCompositeOffsets; ///< Offset of each subdomain in the big composite array 236 unsigned int* SD_SphereCompositeOffsets_SP; ///< like SD_SphereCompositeOffsets, scratch pad (SP) used 237 unsigned int* spheres_in_SD_composite; ///< Big composite array of sphere-subdomain membership 238 239 unsigned int* sphere_owner_SDs; ///< List of owner subdomains for each sphere 240 }; 241 242 // The system is not default-constructible 243 ChSystemGpu_impl() = delete; 244 245 /// Construct Chrono::Gpu system with given sphere radius, density, big domain dimensions and the frame origin. 246 ChSystemGpu_impl(float sphere_rad, float density, float3 boxDims, float3 O); 247 248 /// Create big domain walls out of plane boundary conditions 249 void CreateWallBCs(); 250 251 /// Create an axis-aligned sphere boundary condition 252 size_t CreateBCSphere(float center[3], float radius, bool outward_normal, bool track_forces, float mass); 253 254 /// Create an z-axis aligned cone boundary condition 255 size_t CreateBCConeZ(float cone_tip[3], 256 float slope, 257 float hmax, 258 float hmin, 259 bool outward_normal, 260 bool track_forces); 261 262 /// Create plane boundary condition 263 /// Instead of always push_back, you can select a position in vector to store the BC info: this is for reserved BCs 264 /// (box domain BCs) only. And if the position is SIZE_MAX then the behavior is push_back. 265 size_t CreateBCPlane(float plane_pos[3], float plane_normal[3], bool track_forces, size_t position = SIZE_MAX); 266 267 /// Create customized bc plate 268 size_t CreateCustomizedPlate(float plate_pos_center[3], float plate_normal[3], float hdim_y); 269 270 /// Create an z-axis aligned cylinder boundary condition 271 size_t CreateBCCylinderZ(float center[3], float radius, bool outward_normal, bool track_forces); 272 273 /// Disable a boundary condition by its ID, returns false if the BC does not exist DisableBCbyID(size_t BC_id)274 bool DisableBCbyID(size_t BC_id) { 275 size_t max_id = BC_params_list_SU.size(); 276 if (BC_id >= max_id) { 277 printf("ERROR: Trying to disable invalid BC ID %zu\n", BC_id); 278 return false; 279 } 280 281 if (BC_id <= NUM_RESERVED_BC_IDS - 1) { 282 printf("ERROR: Trying to modify reserved BC ID %zu\n", BC_id); 283 return false; 284 } 285 BC_params_list_UU.at(BC_id).active = false; 286 BC_params_list_SU.at(BC_id).active = false; 287 return true; 288 } 289 290 /// Enable a boundary condition by its ID, returns false if the BC does not exist EnableBCbyID(size_t BC_id)291 bool EnableBCbyID(size_t BC_id) { 292 size_t max_id = BC_params_list_SU.size(); 293 if (BC_id >= max_id) { 294 printf("ERROR: Trying to enable invalid BC ID %zu\n", BC_id); 295 return false; 296 } 297 if (BC_id <= NUM_RESERVED_BC_IDS - 1) { 298 printf("ERROR: Trying to modify reserved BC ID %zu\n", BC_id); 299 return false; 300 } 301 BC_params_list_UU.at(BC_id).active = true; 302 BC_params_list_SU.at(BC_id).active = true; 303 return true; 304 } 305 306 /// Enable a boundary condition by its ID, returns false if the BC does not exist SetBCOffsetFunction(size_t BC_id,const GranPositionFunction & offset_function)307 bool SetBCOffsetFunction(size_t BC_id, const GranPositionFunction& offset_function) { 308 size_t max_id = BC_params_list_SU.size(); 309 if (BC_id >= max_id) { 310 printf("ERROR: Trying to set offset function for invalid BC ID %zu\n", BC_id); 311 return false; 312 } 313 if (BC_id <= NUM_RESERVED_BC_IDS - 1) { 314 printf("ERROR: Trying to modify reserved BC ID %zu\n", BC_id); 315 return false; 316 } 317 BC_offset_function_list.at(BC_id) = offset_function; 318 BC_params_list_UU.at(BC_id).fixed = false; 319 BC_params_list_SU.at(BC_id).fixed = false; 320 return true; 321 } 322 323 /// Prescribe the motion of the big domain, allows wavetank-style simulations setBDWallsMotionFunction(const GranPositionFunction & pos_fn)324 void setBDWallsMotionFunction(const GranPositionFunction& pos_fn) { 325 BDOffsetFunction = pos_fn; 326 BC_offset_function_list.at(BD_WALL_ID_X_BOT) = pos_fn; 327 BC_offset_function_list.at(BD_WALL_ID_X_TOP) = pos_fn; 328 BC_offset_function_list.at(BD_WALL_ID_Y_BOT) = pos_fn; 329 BC_offset_function_list.at(BD_WALL_ID_Y_TOP) = pos_fn; 330 BC_offset_function_list.at(BD_WALL_ID_Z_BOT) = pos_fn; 331 BC_offset_function_list.at(BD_WALL_ID_Z_TOP) = pos_fn; 332 } 333 334 // Copy back the subdomain device data and save it to a file for error checking on the priming kernel 335 //// RADU: is this function going to be implemented?!? 336 ////void checkSDCounts(std::string ofile, bool write_out, bool verbose) const; 337 338 /// Get the max z position of the spheres, allows easier co-simulation. True for getting max Z, false for getting 339 /// minimum Z. 340 double GetMaxParticleZ(bool getMax = true); 341 342 /// Return the total kinetic energy of all particles. 343 float ComputeTotalKE(); 344 345 /// Return the squared sum of the 3 arrays. 346 float computeArray3SquaredSum(std::vector<float, cudallocator<float>>& arrX, 347 std::vector<float, cudallocator<float>>& arrY, 348 std::vector<float, cudallocator<float>>& arrZ, 349 size_t nSpheres); 350 351 /// Return particle position. 352 float3 GetParticlePosition(int nSphere) const; 353 354 /// Set particle position 355 void SetParticlePosition(int nSphere, double3 position); 356 357 /// return absolute velocity 358 float getAbsVelocity(int nSphere); 359 360 // whether or not the particle is fixed 361 bool IsFixed(int nSphere) const; 362 363 /// Return particle linear velocity. 364 float3 GetParticleLinVelocity(int nSphere) const; 365 366 /// Return particle angular velocity. 367 float3 GetParticleAngVelocity(int nSphere) const; 368 369 /// Return particle acceleration 370 float3 GetParticleLinAcc(int nSphere) const; 371 372 /// Return number of particle-particle contacts. 373 int GetNumContacts() const; 374 375 /// Return position of BC plane. 376 float3 GetBCPlanePosition(size_t plane_id) const; 377 378 /// return position of BC sphere 379 float3 GetBCSpherePosition(size_t bc_id) const; 380 381 /// set bc sphere position 382 void SetBCSpherePosition(size_t bc_id, const float3 pos); 383 384 /// set bc sphere vleocity 385 void SetBCSphereVelocity(size_t bc_id, const float3 velo); 386 387 /// return velocity of BC sphere 388 float3 GetBCSphereVelocity(size_t bc_id) const; 389 390 /// Get the reaction forces on a boundary by ID, returns false if the forces are invalid (bad BC ID) 391 bool GetBCReactionForces(size_t BC_id, float3& force) const; 392 393 /// Set initial particle positions. MUST be called only once and MUST be called before initialize 394 void SetParticles(const std::vector<float3>& points, 395 const std::vector<float3>& vels = std::vector<float3>(), 396 const std::vector<float3>& ang_vels = std::vector<float3>()); 397 398 /// set particle velocity, can be called during the simulation 399 void SetParticleVelocity(int id, const double3& velocity); 400 401 void SetBCPlaneRotation(size_t plane_id, double3 rotation_center, double3 rotation_omega); 402 403 /// Advance simulation by duration in user units, return actual duration elapsed 404 /// Requires initialize() to have been called 405 virtual double AdvanceSimulation(float duration); 406 407 /// This method figures out how big a SD is, and how many SDs are going to be necessary 408 /// in order to cover the entire BD. 409 /// Nomenclature: BD: Big domain, SD: Sub-domain. 410 void partitionBD(); 411 412 /// Copy constant sphere data to device 413 void copyConstSphereDataToDevice(); 414 415 /// Reset binning and broadphase info 416 void resetBroadphaseInformation(); 417 /// Reset sphere accelerations 418 void resetSphereAccelerations(); 419 420 /// Reset sphere-wall forces 421 void resetBCForces(); 422 423 /// Collect all the sphere data into the member struct 424 void packSphereDataPointers(); 425 426 /// Run the first sphere broadphase pass to get things started 427 void runSphereBroadphase(); 428 429 /// Wrap the device helper function 430 int3 getSDTripletFromID(unsigned int SD_ID) const; 431 432 /// Create a helper to do sphere initialization 433 void initializeSpheres(); 434 435 /// Sorts particle positions spatially in order to improve memory efficiency 436 void defragment_initial_positions(); 437 438 /// Sorts the user-provided contact history array in the order determined by defragment_initial_positions(), 439 /// if that is called 440 void defragment_friction_history(unsigned int history_offset); 441 442 /// Setup sphere data, initialize local coords 443 void setupSphereDataStructures(); 444 445 /// Helper function to convert a position in UU to its SU representation while also changing data type 446 template <typename T1, typename T2> convertToPosSU(T2 val)447 T1 convertToPosSU(T2 val) { 448 return (T1)(val / gran_params->LENGTH_UNIT); 449 } 450 451 /// convert all BCs from UU to SU 452 void convertBCUnits(); 453 454 /// Max velocity of all particles in system 455 float get_max_vel() const; 456 457 /// Get the maximum stiffness term in the system 458 virtual double get_max_K() const; 459 460 /// This method defines the mass, time, length Simulation Units. It also sets several other constants that enter the 461 /// scaling of various physical quantities set by the user. 462 virtual void switchToSimUnits(); 463 464 /// combine material properties of two types to get effective ones 465 void combineMaterialSurface(); 466 467 /// Set the position function of a boundary condition and account for the offset 468 void setBCOffset(const BC_type&, 469 const BC_params_t<float, float3>& params_UU, 470 BC_params_t<int64_t, int64_t3>& params_SU, 471 double3 offset_UU); 472 473 /// Update positions of each boundary condition using prescribed functions 474 void updateBCPositions(); 475 476 /// Write particle positions, vels and ang vels to a file stream (based on a format) 477 void WriteRawParticles(std::ofstream& ptFile) const; 478 void WriteCsvParticles(std::ofstream& ptFile) const; 479 void WriteChPFParticles(std::ofstream& ptFile) const; 480 #ifdef USE_HDF5 481 void WriteH5Particles(H5::H5File& ptFile) const; 482 #endif 483 484 /// Write contact info file 485 void WriteContactInfoFile(const std::string& outfilename) const; 486 487 /// Get rolling friction torque between body i and j, return 0 if not in contact 488 float3 getRollingFrictionTorque(int i, int j); 489 490 /// get rolling friction v_rot 491 float3 getRollingVrot(int i, int j); 492 493 /// get rolling characterisitc contact time 494 float getRollingCharContactTime(int i, int j); 495 496 /// Get tangential friction force between body i and j, return 0 if not in contact 497 float3 getSlidingFrictionForce(int i, int j); 498 499 /// Get normal friction force between body i and j, return 0 if not in contact 500 float3 getNormalForce(int i, int j); 501 502 /// get list of neighbors in contact with particle ID 503 void getNeighbors(int ID, std::vector<int>& neighborList); 504 505 /// Rough estimate of the total amount of memory used by the system. 506 size_t EstimateMemUsage() const; 507 508 // Conversion factors from SU to UU 509 /// 1 / C_L. Any length expressed in SU is a multiple of SU2UU 510 double LENGTH_SU2UU; 511 /// 1 / C_T. Any time quantity in SU is measured as a positive multiple of TIME_SU2UU 512 double TIME_SU2UU; 513 /// 1 / C_M. Any mass quantity is measured as a positive multiple of MASS_SU2UU. 514 double MASS_SU2UU; 515 /// 1 / C_F. Any force quantity is measured as a multiple of FORCE_SU2UU. 516 double FORCE_SU2UU; 517 /// 1 / C_tau. Any torque quantity is measured as a multiple of TORQUE_SU2UU. 518 double TORQUE_SU2UU; 519 /// 1 / C_v. Any velocity quantity is measured as a multiple of VEL_SU2UU. 520 double VEL_SU2UU; 521 522 // Tuning for non-dimensionalization 523 /// Safety factor on simulation time 524 unsigned int psi_T; 525 526 /// Safety factor on space adim 527 unsigned int psi_L; 528 529 /// Fraction of sphere radius which gives an upper bound on the length unit 530 float psi_R; 531 532 /// Holds the sphere and big-domain-related params in unified memory 533 GranParams* gran_params; 534 535 /// Holds system degrees of freedom 536 SphereData* sphere_data; 537 538 /// Contains information about the status of the granular simulator (solver) 539 ChSolverStateData stateOfSolver_resources; 540 541 /// Allows the code to be very verbose for debugging 542 CHGPU_VERBOSITY verbosity; 543 544 /// If dividing the longest box dimension into INT_MAX pieces gives better resolution than the deformation-based 545 /// scaling, do that. 546 bool use_min_length_unit; 547 548 /// If true, on Initialize(), the order of particles will be re-arranged so those in the same SD locate close to 549 /// each other 550 bool defragment_on_start = true; 551 552 /// Bit flags indicating what fields to write out during WriteParticleFile 553 /// Set with the CHGPU_OUTPUT_FLAGS enum 554 unsigned int output_flags; 555 556 /// How to write the output files? 557 /// Default is CSV 558 CHGPU_OUTPUT_MODE file_write_mode; 559 560 /// Number of discrete elements 561 unsigned int nSpheres; 562 /// Number of subdomains 563 unsigned int nSDs; 564 565 // Use CUDA allocator written by Colin Vanden Heuvel 566 // Could hit system performance if there's not a lot of RAM 567 // Makes somewhat faster memcpys 568 /// Store X positions relative to owner subdomain in unified memory 569 std::vector<int, cudallocator<int>> sphere_local_pos_X; 570 /// Store Y positions relative to owner subdomain in unified memory 571 std::vector<int, cudallocator<int>> sphere_local_pos_Y; 572 /// Store Z positions relative to owner subdomain in unified memory 573 std::vector<int, cudallocator<int>> sphere_local_pos_Z; 574 /// Store X velocity in unified memory 575 std::vector<float, cudallocator<float>> pos_X_dt; 576 /// Store Y velocity in unified memory 577 std::vector<float, cudallocator<float>> pos_Y_dt; 578 /// Store Z velocity in unified memory 579 std::vector<float, cudallocator<float>> pos_Z_dt; 580 581 /// Store X angular velocity (axis and magnitude in one vector) in unified memory 582 std::vector<float, cudallocator<float>> sphere_Omega_X; 583 /// Store Y angular velocity (axis and magnitude in one vector) in unified memory 584 std::vector<float, cudallocator<float>> sphere_Omega_Y; 585 /// Store Z angular velocity (axis and magnitude in one vector) in unified memory 586 std::vector<float, cudallocator<float>> sphere_Omega_Z; 587 588 /// Store X acceleration in unified memory 589 std::vector<float, cudallocator<float>> sphere_acc_X; 590 /// Store Y acceleration in unified memory 591 std::vector<float, cudallocator<float>> sphere_acc_Y; 592 /// Store Z acceleration in unified memory 593 std::vector<float, cudallocator<float>> sphere_acc_Z; 594 595 /// Store X angular acceleration (axis and magnitude in one vector) in unified memory 596 std::vector<float, cudallocator<float>> sphere_ang_acc_X; 597 /// Store Y angular acceleration (axis and magnitude in one vector) in unified memory 598 std::vector<float, cudallocator<float>> sphere_ang_acc_Y; 599 /// Store Z angular acceleration (axis and magnitude in one vector) in unified memory 600 std::vector<float, cudallocator<float>> sphere_ang_acc_Z; 601 602 /// X acceleration history used for the Chung integrator in unified memory 603 std::vector<float, cudallocator<float>> sphere_acc_X_old; 604 /// Y acceleration history used for the Chung integrator in unified memory 605 std::vector<float, cudallocator<float>> sphere_acc_Y_old; 606 /// Z acceleration history used for the Chung integrator in unified memory 607 std::vector<float, cudallocator<float>> sphere_acc_Z_old; 608 /// X angular acceleration history used for the Chung integrator in unified memory 609 std::vector<float, cudallocator<float>> sphere_ang_acc_X_old; 610 /// Y angular acceleration history used for the Chung integrator in unified memory 611 std::vector<float, cudallocator<float>> sphere_ang_acc_Y_old; 612 /// Z angular acceleration history used for the Chung integrator in unified memory 613 std::vector<float, cudallocator<float>> sphere_ang_acc_Z_old; 614 615 /// Fixity of each sphere 616 std::vector<not_stupid_bool, cudallocator<not_stupid_bool>> sphere_fixed; 617 618 /// A buffer array that can store any derived quantity from particles. When users quarry a quantity (such as kinetic 619 /// energy using GetParticleKineticEnergy), this array is used to store that particle-wise quantity, and then 620 /// potentially reduced via CUB. 621 std::vector<float, cudallocator<float>> sphere_stats_buffer; 622 623 /// Set of contact partners for each sphere. Only used in frictional simulations 624 std::vector<unsigned int, cudallocator<unsigned int>> contact_partners_map; 625 /// Whether the frictional contact at an index is active 626 std::vector<not_stupid_bool, cudallocator<not_stupid_bool>> contact_active_map; 627 /// Tracks the tangential history vector for a given contact pair. Only used in multistep friction 628 std::vector<float3, cudallocator<float3>> contact_history_map; 629 /// Tracks the duration of contact between contact pairs. Only used in multistep friction 630 std::vector<float, cudallocator<float>> contact_duration; 631 /// Tracks the normal contact force for a given contact pair 632 std::vector<float3, cudallocator<float3>> normal_contact_force; 633 /// Tracks the tangential contact force for a given contact pair 634 std::vector<float3, cudallocator<float3>> tangential_friction_force; 635 /// Tracks the rolling resistance for a given contact pair 636 std::vector<float3, cudallocator<float3>> rolling_friction_torque; 637 /////////////////////DEBUG PURPOSE/////////////////////////// 638 /// Tracks the characteristic contact time for a given contact pair 639 std::vector<float, cudallocator<float>> char_collision_time; 640 /// Tracks v_rot 641 std::vector<float3, cudallocator<float3>> v_rot_array; 642 643 /// X gravity in user units 644 float X_accGrav; 645 /// Y gravity in user units 646 float Y_accGrav; 647 /// Z gravity in user units 648 float Z_accGrav; 649 650 /// Entry "i" says how many spheres touch subdomain i 651 std::vector<unsigned int, cudallocator<unsigned int>> SD_NumSpheresTouching; 652 /// Entry "i" says where spheres touching ith SD are stored in the big composite array (the offset) 653 std::vector<unsigned int, cudallocator<unsigned int>> SD_SphereCompositeOffsets; 654 /// Scratch area, needed to populate data in the big composite array 655 std::vector<unsigned int, cudallocator<unsigned int>> SD_SphereCompositeOffsets_ScratchPad; 656 /// Array with the IDs of the spheres touching each SD associated with the box; organized SD after SD 657 std::vector<unsigned int, cudallocator<unsigned int>> spheres_in_SD_composite; 658 659 /// List of owner subdomains for each sphere 660 std::vector<unsigned int, cudallocator<unsigned int>> sphere_owner_SDs; 661 662 /// User provided timestep in UU 663 float stepSize_UU; 664 665 /// SU (adimensional) value of time step 666 float stepSize_SU; 667 668 /// how is the system integrated through time? 669 CHGPU_TIME_INTEGRATOR time_integrator; 670 671 /// Total time elapsed since beginning of simulation 672 float elapsedSimTime; 673 674 /// Normal stiffness for sphere-to-sphere contact (Hertzian spring constant) 675 double K_n_s2s_UU; 676 677 /// Normal stiffness for sphere-to-wall contact (Hertzian spring constant) 678 double K_n_s2w_UU; 679 680 /// Normal damping for sphere-to-sphere 681 double Gamma_n_s2s_UU; 682 683 /// Normal damping for sphere-to-wall 684 double Gamma_n_s2w_UU; 685 686 /// Tangential stiffness for sphere-to-sphere (Hertzian spring constant 687 double K_t_s2s_UU; 688 689 /// Tangential stiffness for sphere-to-wall (Hertzian spring constant 690 double K_t_s2w_UU; 691 692 /// Tangential damping for sphere-to-sphere 693 double Gamma_t_s2s_UU; 694 695 /// Tangential damping for sphere-to-wall 696 double Gamma_t_s2w_UU; 697 698 /// material based property: youngs modulus, poisson ratio, cor 699 bool use_mat_based = false; 700 701 double YoungsModulus_sphere_UU; 702 double YoungsModulus_wall_UU; 703 double YoungsModulus_mesh_UU; 704 705 double COR_sphere_UU; 706 double COR_wall_UU; 707 double COR_mesh_UU; 708 709 double PoissonRatio_sphere_UU; 710 double PoissonRatio_wall_UU; 711 double PoissonRatio_mesh_UU; 712 713 /// Rolling friction coefficient for sphere-to-sphere 714 double rolling_coeff_s2s_UU; 715 716 /// Rolling friction coefficient for sphere-to-wall 717 /// Units and use are dependent on the rolling friction model used 718 double rolling_coeff_s2w_UU; 719 720 /// Spinning friction coefficient for sphere-to-sphere 721 double spinning_coeff_s2s_UU; 722 723 /// Spinning friction coefficient for sphere-to-wall 724 /// Units and use are dependent on the spinning friction model used 725 double spinning_coeff_s2w_UU; 726 727 /// Store the ratio of the acceleration due to cohesion vs the acceleration due to gravity, makes simple API 728 float cohesion_over_gravity; 729 730 /// Store the ratio of the acceleration due to adhesion vs the acceleration due to gravity 731 float adhesion_s2w_over_gravity; 732 733 /// The reference point for UU to SU local coordinates 734 int64_t3 BD_rest_frame_SU; 735 736 /// List of generalized boundary conditions that constrain sphere motion 737 std::vector<BC_type, cudallocator<BC_type>> BC_type_list; 738 /// Sim-unit (adimensional) details of boundary conditions 739 std::vector<BC_params_t<int64_t, int64_t3>, cudallocator<BC_params_t<int64_t, int64_t3>>> BC_params_list_SU; 740 /// User-unit (dimensional) details of boundary conditions 741 std::vector<BC_params_t<float, float3>, cudallocator<BC_params_t<float, float3>>> BC_params_list_UU; 742 /// Offset motions functions for boundary conditions -- used for moving walls, wavetanks, etc. 743 std::vector<GranPositionFunction> BC_offset_function_list; 744 745 /// User defined radius of the sphere 746 float sphere_radius_UU; 747 /// User defined density of the sphere 748 float sphere_density_UU; 749 750 /// X-length of the big domain; defines the global X axis located at the CM of the box (as default) 751 float box_size_X; 752 /// Y-length of the big domain; defines the global Y axis located at the CM of the box (as default) 753 float box_size_Y; 754 /// Z-length of the big domain; defines the global Z axis located at the CM of the box (as default) 755 float box_size_Z; 756 757 /// XYZ coordinate of the center of the big box domain in the user-defined frame. Default is (0,0,0). 758 float user_coord_O_X; 759 float user_coord_O_Y; 760 float user_coord_O_Z; 761 762 /// User-provided sphere positions in UU 763 std::vector<float3> user_sphere_positions; 764 765 /// User-provided sphere velocities in UU 766 std::vector<float3> user_sphere_vel; 767 768 /// User-provided sphere angular velocities in UU 769 std::vector<float3> user_sphere_ang_vel; 770 771 /// User-provided sphere fixity as bools 772 std::vector<bool> user_sphere_fixed; 773 774 /// User-provided sphere contact pair information 775 std::vector<unsigned int> user_partner_map; 776 777 /// User-provided sphere (contact pair) friction history information in UU 778 std::vector<float3> user_friction_history; 779 780 /// The offset function for the big domain walls 781 GranPositionFunction BDOffsetFunction; 782 783 /// Allow the user to set the big domain to be fixed, ignoring any given position functions 784 bool BD_is_fixed = true; 785 786 public: 787 // Do two things: make the naming nicer and require a const pointer everywhere 788 789 /// Get handle for the gran params that skips namespacing and enforces const-ness 790 typedef const ChSystemGpu_impl::GranParams* GranParamsPtr; 791 792 /// Get handle for the sphere data that skips namespacing and enforces const-ness 793 typedef const ChSystemGpu_impl::SphereData* GranSphereDataPtr; 794 795 friend class ChSystemGpu; 796 friend class ChSystemGpuMesh; 797 }; 798 799 } // namespace gpu 800 } // namespace chrono 801