Triaxial Compression Test in DEM Engine

[Shahriar Karim Shafin]

Can I do a triaxial compression test in Chrono DEM Engine??  I can easily do uniaxial compression but I cant put any confined pressure in y and z directions when trying triaxial tests in chrono DEM Engine. How can I keep the pressure constant in a direction  at the time of compression??

[Ruochun Zhang]

The short answer is yes, but it's not a one-liner. If you already have your setup to compress the material, then rather than directly setting the wall pressure, think of creating a servo-controlled boundary problem. You can use the tracker class to track the walls that you use for compression, then you can get the normal force using that tracker, and divide that by the area of the wall, you get the pressure. Adjust the wall velocity (or position) based on the difference between the current pressure and your target pressure. This creates a feedback loop that drives the system toward the desired pressure.

Let me know if there are further questions.

Ruochun

[Shahriar Karim Shafin]

Thanks a lot, but the problem is how can I keep that pressure constant after getting my desired pressure (If I put an axial strain rate in the z-direction after that, the lateral pressure will increase again according to you)? I wanted to apply an isotropic confining pressure of some value to the specimen in all directions to consolidate it to equilibrium and then an axial strain rate on the top wall (z-axis) while maintaining constant lateral pressure (x and y axis) which is also called Standard Triaxial Test.

[Ruochun Zhang]

Then again, I still think it's a servo-controlled boundary problem, or say you need a control policy implemented (i.e. you keep the lateral pressure constant by keeping the lateral walls under servo control during the axial loading phase). During shear you prescribe the top wall’s z-velocity (strain rate control) while continuously (like updating it every N steps or so) adjusting the x/y wall positions/velocities based on the tracker-measured lateral pressure values​. This should be exactly like what's happening under the hood for a physical test, and this control policy shouldn't be too hard to implement with C++ or Python. It's not already in DEME, which is a shame... But if you can contribute such a demo to the project, that'd be great. I hope I understand your question correctly.

[Shahriar Karim Shafin]

Got ur point. I will try to do that. I have few questions regarding DEME.

1. Can I use extremely small particle sizes for compression using Hertz-Mindlin model. My specimen size is 50nm currently and its consistently showing me error of very high velocity (exceeding max velocity around 2k) at step 0?

2. When I run a simulation, sometimes its using 1 GPU (getting stuck at some point) or sometimes 2 GPUs (runs smoothly then). Then I had to rerun the simulation . Is there any way to solve it? Particle size count is around 30k.

3. In another compaction problem, I am using particles with E = 200 GPa. I have set the time step accordingly, and the particle size is 1 mm. However, the particles are not settling properly. They keep vibrating(underdamped) even after a long settling period (increasing the settling time or reducing CoR to 0 doesn't help). 

In all the demo examples, I have seen that at most 1 GPa is used for Young’s modulus. How can I stabilize the particle bed when using such a high stiffness material?

Thanks for the support.

[Ruochun Zhang]

Here is my take on your questions:

1. You can. Although I doubt how viable 50nm particles are. You can always scale the length unit if needed, for sure, but with or without doing it, the velocity values will be large compared to particle sizes, necessitating extremely small time step sizes (or you have, for example, particle phasing through one another in one time step). The step size requirement may be quite ridiculous for your needs, like 1e-8 or 1e-9. But another thing is, you need to ensure the system is initialized correctly, that no particle is spawned in the walls (analytical walls may be invisible when you visualize), that is a more typical cause for simulations that dies right at the start.

2. Although this may look like a VRAM issue, from my experience, for DEME, this is more like the kT's automatic bin size adaptation policy is failing. It has a tendency to fail for some "outlier" simulation scenarios, and I hope to resolve that in the next major version. You can probably test if this is the case by adding the following two lines before Initialization(), to see if they help:

DEMSim.SetInitBinNumTarget(int(1e7)); // or some other number like 1e6 or 1e5
DEMSim.DisableAdaptiveBinSize();

If this makes no difference, then we'll have to see the code to figure out what is really happening.

3. Again, using such a high E will necessitate extremely small time step sizes. I don't know what you are using (maybe you derived it using Rayleigh time step?), but it's always good to try if even smaller step sizes stabilize it. It should be noted that even if it does, the size may be too small to be viable. Using reduced E in DEM simulations is common and in many use cases it hardly affects the result.

There's also a chance that the number is too big and the floating-point accuracy starts to kick in. So another good thing to try is to scale the length unit so 200G becomes a smaller number. But this is less likely the cause.

Ruochun

[Shahriar Karim Shafin]

For Question 2, the suggested change did not fix the GPU utilization issue. The simulation remains stuck at Frame 0 during the settling phase and consistently uses only 1 GPU.

However, if I reduce the particle size to 1 mm and scale the other parameters accordingly, the simulation sometimes uses 2 GPUs and runs normally. That said, it is not consistent — sometimes it still uses only 1 GPU and gets stuck, and other times it uses 2 GPUs and runs smoothly.

For the current particle size of 1 micrometer, it is always using 1 GPU and remains stuck at Frame 0 during settling. I am attaching the code snippet here.

int main() {

    DEMSolver DEMSim;

    DEMSim.UseFrictionalHertzianModel();

    DEMSim.SetVerbosity(INFO);

    DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV);

    DEMSim.SetOutputContent(OUTPUT_CONTENT::VEL);

    DEMSim.SetMeshOutputFormat(MESH_FORMAT::VTK);

    DEMSim.SetContactOutputContent(

    CNT_TYPE |

    FORCE |

    CNT_POINT |

    TORQUE |

    CNT_WILDCARD |

    OWNER

);

    DEMSim.SetErrorOutAvgContacts(200);

    //DEMSim.SetInitBinNumTarget(int(1e5)); // 1 million bins is plenty for your setup

    //DEMSim.DisableAdaptiveBinSize();

    auto mat_type_container =

        DEMSim.LoadMaterial({{"E", 1e9}, {"nu", 0.3}, {"CoR", 0.3}, {"mu", 0.80}, {"Crr", 0.10}});

    auto mat_type_particle =

        DEMSim.LoadMaterial({{"E", 1e9}, {"nu", 0.3}, {"CoR", 0.1}, {"mu", 0.40}, {"Crr", 0.04}});

    /* auto mat_type_particle2 =

        DEMSim.LoadMaterial({{"E", 1e10}, {"nu", 0.3}, {"CoR", 0.3}, {"mu", 0.40}, {"Crr", 0.01}});

    auto mat_type_particle3 =

        DEMSim.LoadMaterial({{"E", 1e9}, {"nu", 0.3}, {"CoR", 0.3}, {"mu", 0.40}, {"Crr", 0.04}});

    auto mat_type_particle4 =

        DEMSim.LoadMaterial({{"E", 1e10}, {"nu", 0.3}, {"CoR", 0.3}, {"mu", 0.40}, {"Crr", 0.04}});

    auto mat_type_particle5 =

        DEMSim.LoadMaterial({{"E", 1.7e11}, {"nu", 0.3}, {"CoR", 0.1}, {"mu", 0.40}, {"Crr", 0.04}});

    */

   // Define interaction properties between the container and the particle

    DEMSim.SetMaterialPropertyPair("CoR", mat_type_container, mat_type_particle, 0.1);

    DEMSim.SetMaterialPropertyPair("mu", mat_type_container, mat_type_particle, 0.35);

    DEMSim.SetMaterialPropertyPair("Crr", mat_type_container, mat_type_particle, 0.01);

    /*

    DEMSim.SetMaterialPropertyPair("CoR", mat_type_container, mat_type_particle1, 0.1);

    DEMSim.SetMaterialPropertyPair("mu", mat_type_container, mat_type_particle1, 0.35);

    DEMSim.SetMaterialPropertyPair("Crr", mat_type_container, mat_type_particle1, 0.01);

    */

   // ====================== ADDED RESTART FILE PATHS ======================

    path restart_clump_file;

    path restart_contact_file;

    path restart_state_file;

    bool restart_mode = false;

    double restart_time = 0.0;

   

    path out_dir = current_path();

    out_dir /= "DemoOutput_GravityCompaction_Yttria2"; // Changed folder name

    create_directory(out_dir);  

    restart_clump_file = out_dir / "restart_clumps.csv";

    restart_contact_file = out_dir / "restart_contacts.csv";

    restart_state_file = out_dir / "restart_state.txt";

    restart_mode = exists(restart_clump_file);

    // ===== ADD HERE =====

    if (restart_mode)

        std::cout << "Checkpoint detected — resuming simulation\n";

    else

        std::cout << "Starting fresh simulation\n";

    // ====================

   

    string detected_log = findLatestLogFile();

    cout << "Detected log file: " << (detected_log.empty() ? "(none found)" : detected_log) << "\n";

    // === START ADDED CODE ===

    // 1. Pass paths to globals so the emergency handler can see them

    g_out_dir = out_dir;

    g_log_file = detected_log;

    g_src_file = __FILE__;

    // 2. Register the signals

    signal(SIGINT, signal_callback_handler);  // Catches Ctrl+C

    signal(SIGTERM, signal_callback_handler); // Catches termination requests

    // === END ADDED CODE ===

    // Save metadata at beginning

    saveRunMetadata(out_dir, detected_log, __FILE__);

    float world_size = 0.0015f;

    float container_diameter = 1.1e-4f;

    float terrain_density = 5.03e3;

   

    // --- Log-Normal Parameters ---

    float sphere_rad = 1.0e-6f;          // Mean Radius

    float polydispersity = 0.40;        // Standard deviation as % of mean (0.4 = 40%)

    float compaction_speed_1 = -1.0e-4f;

    float compaction_speed_2 = -1.0e-5f;

    float compaction_speed_3 = 1.0e-5f;  

    float step_size = 1.0e-9f; // Time step size for the simulation

    //float fact_radius = 0.9;

    DEMSim.InstructBoxDomainDimension(world_size, world_size, world_size);

    DEMSim.InstructBoxDomainBoundingBC("none", mat_type_container);

    double bottom = 0;

    double top = 2.75e-4;

    // ======================================================================================

    // ## Setup Compaction Plate ##

    // ======================================================================================

    path mesh_path = "/data/lab/boddeti/DEM-Engine/data/mesh/funnel_left.obj";

    if (!exists(mesh_path)) {

        std::cerr << "ERROR: Mesh file not found at: " << mesh_path << "\n";

        return 1;

    }

    auto walls = DEMSim.AddWavefrontMeshObject(mesh_path.string(), mat_type_container);

    float3 move = make_float3(0.0001375, 0.00, 0 - 2 * sphere_rad);  // z

    float4 rot = make_float4(0.7071, 0, 0, 0.7071);

    walls->Scale(make_float3(0.0022, 0.000275, 0.0055));

    walls->Move(move, rot);

    walls->SetFamily(2);

    DEMSim.SetFamilyFixed(2);

    auto cylinder = DEMSim.AddExternalObject();

    cylinder->AddCylinder(make_float3(0), make_float3(0, 0, 1), 1.8 * container_diameter / 2., mat_type_container, 0);

    cylinder->SetFamily(10);

    DEMSim.SetFamilyFixed(10);

    auto compactor = DEMSim.AddWavefrontMeshObject(mesh_path.string(), mat_type_container);

    move = make_float3(0.0001375, 0.00, top + 2 * sphere_rad);  // z

    rot = make_float4(0.7071, 0, 0, 0.7071);

    compactor->Scale(make_float3(0.0022, 0.000275, 0.0055));

    compactor->Move(move, rot);

    compactor->SetFamily(20);

    DEMSim.SetFamilyFixed(20);

    const std::string vel_z_str_1 = to_string_with_precision(compaction_speed_1);

    DEMSim.SetFamilyPrescribedLinVel(21, "0", "0", vel_z_str_1);                                            

    const std::string vel_z_str_2 = to_string_with_precision(compaction_speed_2);

    DEMSim.SetFamilyPrescribedLinVel(22, "0", "0", vel_z_str_2);

    const std::string vel_z_str_up = to_string_with_precision(compaction_speed_3);

    DEMSim.SetFamilyPrescribedLinVel(23, "0", "0", vel_z_str_up);

    auto compactor_tracker = DEMSim.Track(compactor);

    auto bottom_tracker = DEMSim.Track(walls);

    auto cylinder_tracker = DEMSim.Track(cylinder);

    // ======================================================================================

    // ## Setup Granular Bed (Log-Normal Distribution) ##

    // ======================================================================================

    std::vector<std::shared_ptr<DEMClumpTemplate>> templates_terrain;

    std::vector<double> template_weights; // Stores the probability of selecting each template

    std::vector<double> bin_radii;

   

    // 1. Calculate Log-Normal Parameters from Arithmetic Mean and StdDev

    //    Arithmetic Mean = sphere_rad

    //    Arithmetic StdDev = sphere_rad * polydispersity

    double mean = sphere_rad;

    double stddev = sphere_rad * polydispersity;

   

    // Convert arithmetic statistics to Log-scale statistics (mu_ln, sigma_ln)

    double sigma_ln_sq = std::log(1.0 + (stddev * stddev) / (mean * mean));

    double sigma_ln = std::sqrt(sigma_ln_sq);

    double mu_ln = std::log(mean) - 0.5 * sigma_ln_sq;

    std::cout << "Generating Log-Normal Templates...\n";

    std::cout << "Target Mean Rad: " << mean << ", Target StdDev: " << stddev << "\n";

   

    // 2. Generate Discrete Bins (Templates)

    // We create 30 bins to approximate the continuous curve.

    int num_bins = 15;

    std::vector<int> bin_counts(num_bins, 0); // Counts particles per size

    double min_r_cutoff = mean * 0.5; // Cutoff tails to prevent numerically unstable tiny particles

    double max_r_cutoff = mean * 1.5; // Cutoff tails to prevent huge particles

    double r_step = (max_r_cutoff - min_r_cutoff) / (num_bins - 1);

    double max_generated_radius = 0.0; // We need this for the sampler spacing

    for (int i = 0; i < num_bins; i++) {

        double r = min_r_cutoff + i * r_step;

        bin_radii.push_back(r);

       

        // Calculate Log-Normal PDF value for radius 'r'

        double pdf_val = (1.0 / (r * sigma_ln * std::sqrt(2.0 * math_PI))) * std::exp( -std::pow(std::log(r) - mu_ln, 2) / (2.0 * sigma_ln_sq) );

        // 3. APPLY VOLUME WEIGHTING (The Critical Fix)

        // We multiply by r^3 so that we sample based on mass contribution, not just count.

        // double volume_weight = r * r * r;

       

        // double combined_weight = pdf_val * volume_weight;

        double vol = (4.0 / 3.0) * math_PI * r * r * r;

        double mass = terrain_density * vol;

        auto new_template = DEMSim.LoadSphereType(mass, r, mat_type_particle);

       

        // ADDED FIX: Assign names to templates so they can be read by CSV loader

        char t_name[20];

        sprintf(t_name, "%04d", i);

        new_template->AssignName(std::string(t_name));

       

        templates_terrain.push_back(new_template);

        template_weights.push_back(pdf_val); // Weight by mass contribution

        if (r > max_generated_radius) max_generated_radius = r;

    }

    std::cout << "Loaded " << templates_terrain.size() << " log-normal templates.\n";

    std::cout << "Radius Range: [" << min_r_cutoff << " to " << max_generated_radius << "]\n";

    // ------------------------------------------

    //std::random_device rd;

    std::mt19937 gen(42); // Fixed seed for reproducibility

   

    // Create a Discrete Distribution based on the calculated weights

    // This allows us to pick a template index with the correct statistical probability

    std::discrete_distribution<> dist_lognormal(template_weights.begin(), template_weights.end());

    // ----------------------------

    // (C) HCP Sampler

    // ----------------------------

    float safe_spacing = 2.0f * max_generated_radius;

    HCPSampler sampler(safe_spacing);    

    // ------------------------------------------------

    // (D) BOX sampling region dimensions

    // ------------------------------------------------

    float fill_height = top - bottom;

    float sample_z = bottom + safe_spacing;        

    float halfwidth_x = container_diameter / 2.0f - safe_spacing;

    float halfwidth_y = container_diameter / 2.0f - safe_spacing;

    unsigned int num_particles = 0;

    double total_solid_volume = 0.0;

    std::cout << "Sampling particles in BOX layers...\n";

    if (restart_mode) {

        std::cout << "\n=== RESTART MODE DETECTED ===\n";

        // ADDED FIX: Use DEME built-in loaders matching template names

        auto loaded_xyz = DEMSim.ReadClumpXyzFromCsv(restart_clump_file.string());

        auto loaded_quat = DEMSim.ReadClumpQuatFromCsv(restart_clump_file.string());

        std::vector<float3> in_xyz;

        std::vector<float4> in_quat;

        std::vector<std::shared_ptr<DEMClumpTemplate>> in_types;

        total_solid_volume = 0.0;

        for (int i = 0; i < templates_terrain.size(); i++) {

            char t_name[20];

            sprintf(t_name, "%04d", i);

           

            auto this_type_xyz = loaded_xyz[std::string(t_name)];

            auto this_type_quat = loaded_quat[std::string(t_name)];

           

            for (size_t j = 0; j < this_type_xyz.size(); j++) {

                in_xyz.push_back(this_type_xyz[j]);

                in_quat.push_back(this_type_quat[j]);

                in_types.push_back(templates_terrain[i]);

               

                // Recompute solid volume

                double r = bin_radii[i];

                total_solid_volume += (4.0/3.0) * math_PI * r*r*r;

            }

        }

        DEMClumpBatch restart_batch(in_xyz.size());

        restart_batch.SetTypes(in_types);

        restart_batch.SetPos(in_xyz);

        restart_batch.SetOriQ(in_quat);

        restart_batch.SetFamily(1);

        if (exists(restart_contact_file)) {

            auto pairs = DEMSim.ReadContactPairsFromCsv(restart_contact_file.string());

            auto wcs   = DEMSim.ReadContactWildcardsFromCsv(restart_contact_file.string());

            restart_batch.SetExistingContacts(pairs);

            restart_batch.SetExistingContactWildcards(wcs);

        }

        DEMSim.AddClumps(restart_batch);

        num_particles = in_xyz.size();

        std::cout << "Restart loaded particles: " << num_particles << std::endl;

        std::cout << "Recomputed solid volume: " << total_solid_volume << std::endl;

    }

    else {

        // -------------------------------------------------------

        // (E) Layer-by-layer box sampling

        // -------------------------------------------------------

        while (sample_z < bottom + fill_height) {

            float3 sample_center = make_float3(0, 0, sample_z);

            auto input_xyz = sampler.SampleBox(

                sample_center,

                make_float3(halfwidth_x, halfwidth_y, 0.0f)

            );

            // Assign template based on Hybrid Distribution

            std::vector<std::shared_ptr<DEMClumpTemplate>> template_used(input_xyz.size());

            for (unsigned int i = 0; i < input_xyz.size(); i++) {

                int selected_bin_index = dist_lognormal(gen);

                template_used[i] = templates_terrain[selected_bin_index];

                bin_counts[selected_bin_index]++; // Uncomment if tracking counts

            }

            // Add these particles to DEM

            auto added = DEMSim.AddClumps(template_used, input_xyz);

            added->SetFamily(1);  

            num_particles += input_xyz.size();

            sample_z += safe_spacing;  

        }

    }

    // ======================================================================================

    // ## Post-Generation Calculations & Reports ##

    // ======================================================================================

    if (!restart_mode) {

        total_solid_volume = 0.0;

        for (int i = 0; i < num_bins; i++) {

            if (bin_counts[i] > 0) {

                double r = bin_radii[i];

                double vol_per_particle = (4.0 / 3.0) * math_PI * r * r * r;

                total_solid_volume += bin_counts[i] * vol_per_particle;

            }

        }

    }

    if (!restart_mode) {

        std::cout << "\n============================================\n";

        std::cout << "        PARTICLE SIZE DISTRIBUTION REPORT    \n";

        std::cout << "============================================\n";

        std::cout << "Bin ID | Radius (m)  | Radius (mm) | Count  \n";

        std::cout << "-------|-------------|-------------|--------\n";

        std::ofstream hist_file;

        hist_file.open(out_dir / "particle_size_histogram.csv");

        hist_file << "Radius_m,Count\n";

        for (int i = 0; i < num_bins; i++) {

            if (bin_counts[i] > 0) {

                std::cout << std::setw(6) << i << " | "

                          << std::scientific << std::setprecision(3) << bin_radii[i] << " | "

                          << std::fixed << std::setprecision(4) << bin_radii[i] * 1000.0 << "    | "

                          << bin_counts[i] << "\n";

                hist_file << bin_radii[i] << "," << bin_counts[i] << "\n";

            }

        }

        hist_file.close();

        std::cout << "============================================\n";

        std::cout << "Total num of particles: " << num_particles << std::endl;

    }

    auto max_z_finder = DEMSim.CreateInspector("clump_max_z");

    auto min_z_finder = DEMSim.CreateInspector("clump_min_z");

    auto max_v_finder = DEMSim.CreateInspector("clump_max_absv");

    // <-- ADDED: Set the extra margin for particle family 1

   // DEMSim.SetFamilyExtraMargin(1, fact_radius * sphere_rad);

    DEMSim.SetInitTimeStep(step_size);

    DEMSim.SetMaxVelocity(30.);

    DEMSim.SetGravitationalAcceleration(make_float3(0, 0, -9.81));

    //DEMSim.SetErrorOutVelocity(100.);

    DEMSim.Initialize(true); // <-- MODIFIED: Do a dry run

    // <-- ADDED: Setup wildcards and stress-strain output file

    std::cout << "Initial number of contacts: " << DEMSim.GetNumContacts() << std::endl;

 //   DEMSim.SetFamilyContactWildcardValueBoth(1, "initialLength", 0.0);

 //   DEMSim.SetFamilyContactWildcardValueBoth(1, "unbroken", 0.0);

    path stress_file_path = out_dir / "stress_strain_data.csv";

    std::ofstream stress_strain_file;

    if (restart_mode) {

        // Continue writing

        stress_strain_file.open(stress_file_path, std::ios_base::app);

    }

    else {

        // Fresh run → overwrite old file

        stress_strain_file.open(stress_file_path, std::ios_base::out);

        stress_strain_file << "Strain,Stress_Pa,BottomStress_Pa,radial_stress_Pa,Porosity,Void_Ratio,Coordination_Number\n";

    }

    // <-- END ADDED BLOCK

    // =================================================================================================

    // ## Settle the Bed ##

    // =================================================================================================

    float settle_time = 1.0;

    float sim_time = 2.5; // Total simulation time

    unsigned int fps = 200;

    float frame_time = 1.0 / fps;

    unsigned int currframe = 0;

    double stress = 0;

    bool status_1 = true;

    bool status_2 = true;

    bool reversed_back = false;

    double L0 = 0.0;   // <-- declare BEFORE this block

    double restart_plate_z = top + 2 * sphere_rad; // ADDED FIX: Plate Z position tracker

    if (restart_mode && exists(restart_state_file)) {

        std::ifstream state_in(restart_state_file);

        std::string key;

        while (state_in >> key) {

            if (key == "frame") state_in >> currframe;

            else if (key == "time") state_in >> restart_time;

            else if (key == "status1") state_in >> status_1;

            else if (key == "status2") state_in >> status_2;

            else if (key == "reversed") state_in >> reversed_back;

            else if (key == "L0") state_in >> L0;  

            else if (key == "plate_z") state_in >> restart_plate_z; // ADDED FIX

        }

        std::cout << "Restarting from frame " << currframe

                << " at time " << restart_time

                << " with L0 = " << L0 << std::endl;

    }

    if (!restart_mode) {

        std::cout << "Settling bed for " << settle_time << " seconds..." << std::endl;

        for (float t = 0; t < settle_time; t += frame_time) {

            float max_v = max_v_finder->GetValue();

            std::cout << "Frame: " << currframe

                      << " | Time: " << t

                      << " | Max Vel: " << max_v << " m/s"

                      << std::endl;

            char filename[100], meshfilename[100];

            char cnt_filename[100];

            sprintf(filename, "DEMdemo_output_%04d.csv", currframe);

            sprintf(meshfilename, "DEMdemo_mesh_%04d.vtk", currframe);

            sprintf(cnt_filename, "DEMdemo_contact_%04d.csv", currframe);

            DEMSim.WriteSphereFile(out_dir / filename);

            DEMSim.WriteMeshFile(out_dir / meshfilename);

            DEMSim.WriteContactFile(out_dir / cnt_filename);

            currframe++;

            DEMSim.DoDynamicsThenSync(frame_time);

        }

    }

   

    double settled_terrain_max_z = max_z_finder->GetValue();

    double settled_terrain_min_z = min_z_finder->GetValue();

//  float max_v = max_v_finder->GetValue();

    std::cout << "Settled terrain height: " << settled_terrain_max_z << std::endl;

    // <-- ADDED: Get L0 and activate bonds

    if (!restart_mode) {

        L0 = settled_terrain_max_z - settled_terrain_min_z + 2 * sphere_rad;

        std::cout << "Maximum Bed Height: " << L0 << std::endl;

//      std::cout << "Max velocity of any point in simulation is " << max_v << std::endl;

        float dz = settled_terrain_max_z - top;

        compactor_tracker->SetPos(make_float3(0.00, 0.00, dz));

    }

    else {

        // Restore plate position based on saved state

        compactor_tracker->SetPos(make_float3(0.00, 0.00, restart_plate_z)); // ADDED FIX

       

        // Robust Kinematic State Restoration

        if (reversed_back) {

            // Plate was moving upward

            DEMSim.ChangeFamily(20, 23);

            DEMSim.ChangeFamily(21, 23);

            DEMSim.ChangeFamily(22, 23);

        }

        else if (!status_2) {

            // Speed 2 was active

            DEMSim.ChangeFamily(20, 22);

            DEMSim.ChangeFamily(21, 22);

        }

        else if (!status_1) {

            // Speed 1 was active

            DEMSim.ChangeFamily(20, 21);

        }

        // If status_1 is STILL true, it means compaction hadn't started yet

        // at the time of the checkpoint. It safely remains Family 20.

    }

    /*

    DEMSim.SetFamilyClumpMaterial(1, mat_type_particle2); // <-- MODIFIED: Change material of particles to stiffer one before compaction

    DEMSim.DoDynamicsThenSync(0.3);

    std::cout << "Changing particle stiffness(2) to prepare for compaction..." << std::endl;

 

    DEMSim.SetFamilyClumpMaterial(1, mat_type_particle3); // <-- MODIFIED: Change material of particles to stiffer one before compaction

    DEMSim.DoDynamicsThenSync(0.1);

    std::cout << "Changing particle stiffness(3) to prepare for compaction..." << std::endl;

    DEMSim.SetFamilyClumpMaterial(1, mat_type_particle4); // <-- MODIFIED: Change material of particles to stiffer one before compaction

    DEMSim.DoDynamicsThenSync(0.1);

    std::cout << "Changing particle stiffness(4) to prepare for compaction..." << std::endl;

    DEMSim.SetFamilyClumpMaterial(1, mat_type_particle5); // <-- MODIFIED: Change material of particles to stiffer one before compaction

    DEMSim.DoDynamicsThenSync(0.3);

    std::cout << "Changing particle stiffness(5) to prepare for compaction..." << std::endl;

    */

    // ======================================================================================

    // ## Constant Velocity Compaction ##

    // ======================================================================================

    std::cout << "\n=== STARTING CONSTANT VELOCITY COMPACTION ===" << std::endl;

................. After that compaction phase but its stuck at Frame-0 in the settling phase

[Ruochun Zhang]

That inconsistency part is a bit surprising, because right now DEME does not decide what hardware to use based on the problem; it just uses all available GPUs (which, I would say, is a problem and I'd probably have a new version that allows the user to choose how many and which GPUs to use). It almost sounds like one of the devices can go offline in some cases. So my question is, how do you know how many GPUs it is using? From DEME's "active device" output or nvidia-smi readings?

Ruochun

[Shahriar Karim Shafin]

Got your point. I have given you the screenshot. It's from Nvidia-smi readings

[Ruochun Zhang]

Then I assume that DEME still reports 2 devices, yet one of them shows a low utility rate. That's because one of the threads is simply waiting for the update from the other, which never arrives or arrives slowly. More likely, it's kT idling. You can set the verbosity to DEBUG to see how the two threads are collaborating and whether it is idling at all. I guess they are still working but kT is extremely slow for some reason.

You mentioned that this happens when the particle size is scaled down, but did you scale down the world size with it? The kT's automatic bin adaptation generally works well but if the world size is, say, 3 orders of magnitude larger than the active simulation region, it may have some difficulty. If that is not the cause, then it's not easy to guess without having a minimal reproducible example.

Ruochun

[Shahriar Karim Shafin]

Thanks, you are exactly right. kT is extremely slow( runtime around 220s per step). That's why another gpu (dT) is idle. My current particle size is 1 micrometer and timestep 1ns. Particle no is around 130K. I can't change the particle properties (radius, Young's Moduli etc).  Is there any way to make it faster??

[Ruochun Zhang]

Again, to resolve that, we need to have a concrete idea of why kT is so slow. For moderately sized simulations like you have, most typically, it's that for some reason, the simulation entities are moving at high velocity, so the solver considers all n² potential contacts possible or something like that (in which case the physics went bad anyway, regardless of numerical implementation). There are two metrics we should inspect: 1. What is the maximum system velocity? 2. What is the solver-reported system memory usage?

[Shahriar Karim Shafin]

Thanks. I have solved the problem.  The code is working well now, even for nanometer-sized particles, and the kT runtime has also decreased significantly. The earlier issue was caused by precision problems when using SI units as the base system, since the mass and moment of inertia were becoming almost zero (around 20 decimal precision). I converted the entire system to micrometer–microgram–microsecond units, and it is working properly now.

I really appreciate your support.

[Ruochun Zhang]

Then there should be a warning message indicating the mass-related issue. I did mention this being one of the possibilities in the earlier message, but I have to admit, it's been "every time we suspected it but it never turned out to be the true cause" type of problem. It's the first time it actually made a difference, perhaps the mass and MOI values you used finally reached the threshold. In any case, it's good that you fixed it and we can rest easy.

Thank you,

Ruochun