Merge pull request #45 from uwsbel/Mesh_Particles

Mesh particles performance update

Author: Ruochun

cmake/CudaSupportedArchitectures.cmake

@@ -31,7 +31,8 @@ function(cuda_supported_architectures)
 	set(cu10 30 35 50 52 60 61 70 72 75)
 	set(cu11 35 50 52 60 61 70 72 75 80)
 	set(cu11_x 35 50 52 60 61 70 72 75 80 86)
-	set(cu12_x 50 52 60 61 70 72 75 80 86)
+	set(cu12_x 50 52 60 61 70 72 75 80 86 89 120)
+	set(cu13_x 75 80 86 89 90 100 120 121)
 
 	if (CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER_EQUAL 7)
 		set(CUDASUP_ARCHITECTURES ${cu7} CACHE INTERNAL "")
@@ -60,6 +61,10 @@ function(cuda_supported_architectures)
 	if (CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER_EQUAL 12)
 		set(CUDASUP_ARCHITECTURES ${cu12_x} CACHE INTERNAL "")
 	endif()
+	
+    if (CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER_EQUAL 13)
+        set(CUDASUP_ARCHITECTURES ${cu13_x} CACHE INTERNAL "")
+    endif()
 
 	if (NOT DEFINED CUDASUP_ARCHITECTURES)
 		message(SEND_ERROR "[CUDASUP] Could not determine device architectures supported by the CUDA toolkit!")

src/DEM/API.h

@@ -108,7 +108,7 @@ class DEMSolver {
     /// Set the simulation time manually.
     void SetSimTime(double time);
     /// @brief Set the strategy for auto-adapting time step size (NOT implemented, no effect yet).
-    /// @param type "none" or "max_vel" or "int_diff".
+    /// @param type "none" or "hertz_const" or "max_vel" or "int_diff".
     void SetAdaptiveTimeStepType(const std::string& type);
 
     /// @brief Set the time integrator for this simulator.
@@ -222,6 +222,10 @@ class DEMSolver {
     /// thinckness of the contact `safety' margin. This need not to be large unless the simulation velocity can increase
     /// significantly in one kT update cycle.
     void SetExpandSafetyAdder(float vel) { m_expand_base_vel = vel; }
+    /// @brief Control whether angular velocity contributes to the contact detection margin.
+    /// @details Default is auto: false for pure single-sphere clumps; true when multi-sphere clumps or meshes exist.
+    ///          If disabled, sphere--sphere contact uses only linear velocity when computing damping/friction terms.
+    void SetUseAngularVelocityMargin(bool use);
 
     /// @brief Manually set the current triangle--triangle penetration value in dT.
     /// @details This allows the user to directly control the maxTriTriPenetration value which will ONLY be used in the
@@ -325,7 +329,7 @@ class DEMSolver {
     }
     /// @brief Get the current bin (for contact detection) size. Must be called from synchronized stance.
     /// @return Bin size.
-    double GetBinSize() { return kT->simParams->binSize; }
+    double GetBinSize() { return kT->simParams->dyn.binSize; }
     // NOTE: No need to get binSize from the device, as binSize is only changed on the host
 
     /// @brief Get the current number of bins (for contact detection). Must be called from synchronized stance.
@@ -346,6 +350,19 @@ class DEMSolver {
     /// @brief Set the number of past kT updates that dT will use to calibrate the max future drift limit.
     /// @param n Number of kT updates. Suggest using default.
     void SetCDNumStepsMaxDriftHistorySize(unsigned int n);
+    /// @brief Inflate observed kT lag (in units of dT steps) by this factor, for drift scheduling/calibration only.
+    /// @details Values < 1 can be unsafe. Default is 1.1.
+    void SetCDFutureDriftEffDriftSafetyFactor(float factor) {
+        future_drift_eff_drift_safety_factor = clampBetween(factor, 1.0f, 10.0f);
+    }
+    /// @brief Set bounds (as a fraction of max drift) that clamp dT's waiting behavior before sending a new kT work order.
+    /// @details Default is lower=0, upper=1 (send as late as safely possible). `upper=0` makes dT send immediately
+    /// (close to the old behavior, but still with the safety factor). Values must satisfy 0 <= lower <= upper <= 1.
+    void SetCDFutureDriftSendBounds(float lower_ratio, float upper_ratio) {
+        future_drift_send_upper_bound_ratio = clampBetween(upper_ratio, 0.0f, 1.0f);
+        future_drift_send_lower_bound_ratio =
+            clampBetween(lower_ratio, 0.0f, static_cast<float>(future_drift_send_upper_bound_ratio));
+    }
     /// @brief Get the current update frequency used by the solver.
     /// @return The current update frequency.
     float GetUpdateFreq() const;
@@ -1376,6 +1393,13 @@ class DEMSolver {
     /// Show the wall time and percentages of wall time spend on various solver tasks.
     void ShowTimingStats();
 
+    /// Enable/disable GPU event-based timing (StartGpuTimer/StopGpuTimer/AccumulateGpuTimer). Disabling this can
+    /// improve performance by avoiding per-step cudaEventRecord/cudaEventElapsedTime overhead. Call only when the
+    /// solver is not actively running (best: right after constructing the solver, or after DoDynamicsThenSync).
+    void SetGPUTimersEnabled(bool enabled);
+    /// Query whether GPU event-based timing is enabled.
+    bool GetGPUTimersEnabled() const { return m_gpu_timers_enabled; }
+
     /// Reset the collaboration stats between dT and kT back to the initial value (0). You should call this if you want
     /// to start over and re-inspect the stats of the new run; otherwise, it is generally not needed, you can go ahead
     /// and destroy DEMSolver.
@@ -1471,28 +1495,28 @@ class DEMSolver {
     void PrintKinematicScratchSpaceUsage() const { kT->printScratchSpaceUsage(); }
 
     /// Let dT do this call and return the reduce value of the inspected quantity.
-    float dTInspectReduce(const std::shared_ptr<jitify::Program>& inspection_kernel,
+    float dTInspectReduce(const std::shared_ptr<JitHelper::CachedProgram>& inspection_kernel,
                           const std::string& kernel_name,
                           INSPECT_ENTITY_TYPE thing_to_insp,
                           CUB_REDUCE_FLAVOR reduce_flavor,
                           bool all_domain,
                           DualArray<scratch_t>& reduceResArr,
                           DualArray<scratch_t>& reduceRes);
-    float* dTInspectNoReduce(const std::shared_ptr<jitify::Program>& inspection_kernel,
+    float* dTInspectNoReduce(const std::shared_ptr<JitHelper::CachedProgram>& inspection_kernel,
                              const std::string& kernel_name,
                              INSPECT_ENTITY_TYPE thing_to_insp,
                              CUB_REDUCE_FLAVOR reduce_flavor,
                              bool all_domain,
                              DualArray<scratch_t>& reduceResArr,
                              DualArray<scratch_t>& reduceRes);
-    float dTInspectReduceDevice(const std::shared_ptr<jitify::Program>& inspection_kernel,
+    float dTInspectReduceDevice(const std::shared_ptr<JitHelper::CachedProgram>& inspection_kernel,
                                 const std::string& kernel_name,
                                 INSPECT_ENTITY_TYPE thing_to_insp,
                                 CUB_REDUCE_FLAVOR reduce_flavor,
                                 bool all_domain,
                                 DualArray<scratch_t>& reduceResArr,
                                 DualArray<scratch_t>& reduceRes);
-    float* dTInspectNoReduceDevice(const std::shared_ptr<jitify::Program>& inspection_kernel,
+    float* dTInspectNoReduceDevice(const std::shared_ptr<JitHelper::CachedProgram>& inspection_kernel,
                                    const std::string& kernel_name,
                                    INSPECT_ENTITY_TYPE thing_to_insp,
                                    CUB_REDUCE_FLAVOR reduce_flavor,
@@ -1599,6 +1623,9 @@ class DEMSolver {
     // The `base' velocity we always consider entities to have, when determining the thickness of the margin to add for
     // contact detection.
     float m_expand_base_vel = 3.f;
+    // Whether angular velocity contributes to the contact margin (auto-detected if not user-set).
+    bool m_use_angvel_margin = true;
+    bool m_use_angvel_margin_user_set = false;
 
     // The method of determining the thickness of the margin added to CD
     // Default is using a max_vel inspector of the clumps to decide it
@@ -1664,12 +1691,18 @@ class DEMSolver {
     float max_drift_ahead_of_avg_drift = 4.;
     float max_drift_multiple_of_avg_drift = 1.05;
     unsigned int max_drift_gauge_history_size = 200;
+    float future_drift_eff_drift_safety_factor = 1.1;
+    float future_drift_send_upper_bound_ratio = 1.0;
+    float future_drift_send_lower_bound_ratio = 0.0;
 
     // See SetNoForceRecord
     bool no_recording_contact_forces = false;
     // See SetCollectAccRightAfterForceCalc
     bool collect_force_in_force_kernel = false;
 
+    // Controls event-based GPU timing for kT/dT timers.
+    bool m_gpu_timers_enabled = true;
+
     // Error-out avg num contacts
     float threshold_error_out_num_cnts = 300.;
 

src/DEM/APIPrivate.cpp

@@ -295,7 +295,7 @@ void DEMSolver::jitifyKernels() {
     equipIntegrationScheme(m_subs);
     equipKernelIncludes(m_subs);
 
-    // Jitify may require a defined device to derive the arch
+        // Jitify may require a defined device to derive the arch
     std::thread kT_build([&]() {
         DEME_GPU_CALL(cudaSetDevice(kT->streamInfo.device));
         kT->jitifyKernels(m_subs, m_jitify_options);
@@ -304,12 +304,6 @@ void DEMSolver::jitifyKernels() {
     std::thread dT_build([&]() {
         DEME_GPU_CALL(cudaSetDevice(dT->streamInfo.device));
         dT->jitifyKernels(m_subs, m_jitify_options);
-
-        // Now, inspectors need to be jitified too... but the current design jitify inspector kernels at the first time
-        // they are used. for (auto& insp : m_inspectors) {
-        //     insp->Initialize(m_subs);
-        // }
-
         // Solver system's own max vel inspector should be init-ed. Don't bother init-ing it while using, because it is
         // called at high frequency, let's save an if check. Forced initialization (since doing it before system
         // completes init).
@@ -319,8 +313,20 @@ void DEMSolver::jitifyKernels() {
         m_approx_angvel_func->Initialize(m_subs, m_jitify_options, true);
         dT->approxAngVelFunc = m_approx_angvel_func;
     });
+
     kT_build.join();
     dT_build.join();
+
+    // Eagerly initialize user-created inspectors so their kernels compile before first use
+    for (auto& insp : m_inspectors) {
+        if (insp) {
+            insp->Initialize(m_subs, m_jitify_options, true);
+            if (insp->inspection_kernel) {
+                insp->inspection_kernel->kernel(insp->kernel_name).instantiate();
+            }
+        }
+    }
+
 }
 
 void DEMSolver::getContacts_impl(std::vector<bodyID_t>& idA,
@@ -1145,7 +1151,13 @@ void DEMSolver::setSolverParams() {
     dT->solverFlags.upperBoundFutureDrift = upper_bound_future_drift;
     dT->solverFlags.targetDriftMoreThanAvg = max_drift_ahead_of_avg_drift;
     dT->solverFlags.targetDriftMultipleOfAvg = max_drift_multiple_of_avg_drift;
-    dT->accumStepUpdater.SetCacheSize(max_drift_gauge_history_size);
+    kT->solverFlags.futureDriftEffDriftSafetyFactor = future_drift_eff_drift_safety_factor;
+    dT->solverFlags.futureDriftEffDriftSafetyFactor = future_drift_eff_drift_safety_factor;
+    kT->solverFlags.futureDriftSendUpperBoundRatio = future_drift_send_upper_bound_ratio;
+    dT->solverFlags.futureDriftSendUpperBoundRatio = future_drift_send_upper_bound_ratio;
+    kT->solverFlags.futureDriftSendLowerBoundRatio = future_drift_send_lower_bound_ratio;
+    dT->solverFlags.futureDriftSendLowerBoundRatio = future_drift_send_lower_bound_ratio;
+    // accumStepUpdater removed; FutureDriftRegulator manages its own window size.
 }
 
 void DEMSolver::setSimParams() {
@@ -1179,20 +1191,71 @@ void DEMSolver::setSimParams() {
             DEME_MAX_WILDCARD_NUM);
     }
     DEME_DEBUG_PRINTF("%u contact wildcards are in the force model.", nContactWildcards);
-
     // Error-out velocity should be no smaller than the max velocity we can expect
     if (threshold_error_out_vel < m_approx_max_vel) {
         // Silently bring down m_approx_max_vel
         m_approx_max_vel = threshold_error_out_vel;
     }
-
+    { // Adaptive Timestep -- currently only Hertz const. 
+        if (adapt_ts_type == ADAPT_TS_TYPE::HERTZ_CONST) {
+            auto sqr = [](double x) { return x * x; };
+            auto effective_E = [&](double E1, double nu1, double E2, double nu2) -> double {
+                if (E1 <= 0.0 || E2 <= 0.0) return 0.0;
+                return 1.0 / ( ((1.0 - sqr(nu1)) / E1) + ((1.0 - sqr(nu2)) / E2) );
+            };  // max effektive E* over all material contacts
+            double Eeff_max = 0.0;
+            for (size_t i = 0; i < m_loaded_materials.size(); ++i) {
+                const auto& A = m_loaded_materials[i]->mat_prop;
+                const double E1  = (A.count("E")  ? (double)A.at("E")  : 0.0);
+                const double nu1 = (A.count("nu") ? (double)A.at("nu") : 0.3);
+                for (size_t j = i; j < m_loaded_materials.size(); ++j) {
+                    const auto& B = m_loaded_materials[j]->mat_prop;
+                    const double E2  = (B.count("E")  ? (double)B.at("E")  : 0.0);
+                    const double nu2 = (B.count("nu") ? (double)B.at("nu") : 0.3);
+                    const double Ee = effective_E(E1, nu1, E2, nu2);
+                    if (Ee > Eeff_max) Eeff_max = Ee;
+                }
+            } // lowest mass
+            double min_mass = std::numeric_limits<double>::infinity();
+            for (double m : m_template_clump_mass)
+                if (m > 0.0 && m < min_mass) min_mass = m;
+            const double r_min = (double)m_smallest_radius;
+            if (std::isfinite(min_mass) && r_min > 0.0 && Eeff_max > 0.0) {
+                const double R_eff = 0.5 * r_min;
+                const double m_eff = 0.5 * min_mass;
+                const double KH = FOUR_OVER_THREE * std::sqrt(0.1) * Eeff_max * R_eff;
+                const double dt_hertz = (PI / (2.0 * N_DT)) * std::sqrt(m_eff / KH);
+                if (dt_hertz > 0.0 && std::isfinite(dt_hertz)) {
+                    m_ts_size = dt_hertz;  // <- set const. timestep
+                    DEME_INFO("Adaptive time step 'hertz_const': dt = %.9g  (m_min=%.6g, r_min=%.6g, E*=%.6g, N_DT=%.1f)",
+                            m_ts_size, min_mass, r_min, Eeff_max, N_DT);
+                } else {
+                    DEME_WARNING("hertz_const erzeugte ungueltigen dt; behalte bisherigen Wert %.7g.", m_ts_size);
+                }
+            } else {
+                DEME_WARNING("hertz_const: fehlende/ungueltige Daten (m_min=%g, r_min=%g, E*=%g); dt bleibt %.7g.",
+                            min_mass, r_min, Eeff_max, m_ts_size);
+            }
+        }
+    }
+    if (!m_use_angvel_margin_user_set) {
+        bool has_multi_sphere_clump = false;
+        for (const auto& radii : m_template_sp_radii) {
+            if (radii.size() > 1) {
+                has_multi_sphere_clump = true;
+                break;
+            }
+        }
+        const bool has_mesh = nTriGM > 0;
+        m_use_angvel_margin = has_multi_sphere_clump || has_mesh;
+    }
     dT->setSimParams(nvXp2, nvYp2, nvZp2, l, m_voxelSize, m_binSize, nbX, nbY, nbZ, m_boxLBF, m_user_box_min,
                      m_user_box_max, G, m_ts_size, m_expand_factor, m_approx_max_vel, m_max_tritri_penetration,
-                     m_expand_safety_multi, m_expand_base_vel, m_force_model->m_contact_wildcards,
+                     m_expand_safety_multi, m_expand_base_vel, m_use_angvel_margin, m_force_model->m_contact_wildcards,
                      m_force_model->m_owner_wildcards, m_force_model->m_geo_wildcards);
     kT->setSimParams(nvXp2, nvYp2, nvZp2, l, m_voxelSize, m_binSize, nbX, nbY, nbZ, m_boxLBF, m_user_box_min,
                      m_user_box_max, G, m_ts_size, m_expand_factor, m_approx_max_vel, m_max_tritri_penetration,
-                     m_expand_safety_multi, m_expand_base_vel, m_force_model->m_contact_wildcards,
+                     m_expand_safety_multi, m_expand_base_vel, m_use_angvel_margin, m_force_model->m_contact_wildcards,
                      m_force_model->m_owner_wildcards, m_force_model->m_geo_wildcards);
 }
 
@@ -1327,7 +1390,10 @@ void DEMSolver::packDataPointers() {
 }
 
 void DEMSolver::migrateSimParamsToDevice() {
+    DEME_GPU_CALL(cudaSetDevice(dT->streamInfo.device));
     dT->simParams.toDevice();
+
+    DEME_GPU_CALL(cudaSetDevice(kT->streamInfo.device));
     kT->simParams.toDevice();
 }
 
@@ -1562,6 +1628,14 @@ inline void DEMSolver::equipForceModel(std::unordered_map<std::string, std::stri
     if (!no_recording_contact_forces) {
         contact_info_write_strat = FORCE_INFO_WRITE_BACK_STRAT();
     }
+    std::string contact_info_clear_strat = " ";
+    if (!no_recording_contact_forces) {
+        contact_info_clear_strat =
+            "granData->contactPointGeometryA[myContactID] = make_float3(0,0,0);"
+            "granData->contactPointGeometryB[myContactID] = make_float3(0,0,0);"
+            "granData->contactForces[myContactID] = make_float3(0,0,0);"
+            "granData->contactTorque_convToForce[myContactID] = make_float3(0,0,0)";
+    }
 
     if (ensure_kernel_line_num) {
         model = compact_code(model);
@@ -1570,6 +1644,7 @@ inline void DEMSolver::equipForceModel(std::unordered_map<std::string, std::stri
         ingredient_acquisition_B = compact_code(ingredient_acquisition_B);
         whether_reduce_in_kernel = compact_code(whether_reduce_in_kernel);
         contact_info_write_strat = compact_code(contact_info_write_strat);
+        contact_info_clear_strat = compact_code(contact_info_clear_strat);
     }
     strMap["_DEMForceModel_;"] = model;
     strMap["_forceModelPrerequisites_;"] = model_prerequisites;
@@ -1590,8 +1665,9 @@ inline void DEMSolver::equipForceModel(std::unordered_map<std::string, std::stri
     strMap["_forceModelContactWildcardWrite_;"] = cnt_wildcard_write_back;
     strMap["_forceModelContactWildcardDestroy_;"] = cnt_wildcard_destroy_record;
 
-    strMap["_forceCollectInPlaceStrat_;"] = whether_reduce_in_kernel;
-    strMap["_contactInfoWrite_;"] = contact_info_write_strat;
+    strMap["_forceCollectInPlaceStrat_"] = whether_reduce_in_kernel;
+    strMap["_contactInfoWrite_"] = contact_info_write_strat;
+    strMap["_contactInfoClear_"] = contact_info_clear_strat;
 
     DEME_DEBUG_PRINTF("Model ingredient definition:\n%s", ingredient_definition.c_str());