Dense packing for polydisperse cases

[Shahriar Karim Shafin]

I want to create a densely packed particle bed with the maximum possible packing density. For the monodisperse case, I can achieve a packing density of approximately 0.74 using HCPSampler. However, for polydisperse cases, I am having difficulty achieving a similarly dense initial packing. If I use HCPSampler with a spacing of 2rmin​, the larger particles overlap, so I have to use 2rmax​ to avoid overlaps. However, this creates significant gaps between particles and results in a much lower packing density. Is there a way to generate an initially jammed or highly dense packing for a polydisperse particle system?

[Dan Negrut]

Hi Shahriar,

 

The short answer is that no sampler setting will get you there: HCPSampler and GridSampler place points on a lattice, and a lattice has exactly one characteristic spacing. There is no lattice that densely packs a polydisperse assembly. So you are not missing an option. You are asking the sampler to do a job it structurally cannot do.

 

Density in a polydisperse bed comes from dynamics, not from sampling. The standard recipe is to sample LOOSE on purpose, exactly as you did with 2*rmax, and then let gravity close the gaps:

 

1. Sample a batch with spacing based on the largest radius, so nothing overlaps at insertion.

2. Drop the batch in and run a couple of seconds of settling.

3. Repeat, adding batches on top of the settled material until the bed is as deep as you need.

4. Save the settled state and reuse it, so you pay this cost once.

 

Those gaps you are worried about are expected and they close during step 2. The sampler's only job is to place non-overlapping seeds.

 

If you are using DEM-Engine, this exact workflow is already written up as a demo series: DEMdemo_GRCPrep_Part1, Part2 and Part3 in src/demo/. It builds a bed of GRC lunar simulant, which is polydisperse in both size and shape, so it is your problem exactly. Part1 uses an HCPSampler at 2.2 times the largest scale (deliberately loose), adds particles in batches with a small downward velocity so they collapse rather than hang, runs about two seconds of settling per batch, and writes a checkpoint at the end. Part2 then replicates that settled bed to build thickness. Copy that pattern rather than reinventing it.

 

On the Chrono side the same idea appears in demo_DEM_repose.cpp and in the vehicle co-simulation granular terrain nodes, which generate in layers, advancing the insertion height by roughly one particle diameter per layer.

 

Two levers if you still want more density after that. First, run the settling phase with friction set very low or to zero and restore the real value afterwards. Friction is what locks a granular assembly into a loose state, so a near-frictionless settle reaches a noticeably denser packing. Second, finish with a light compression plate or a few taps rather than relying on gravity alone.

 

One expectation worth resetting, and I think this is the real answer to your question. The 0.74 you are getting in the monodisperse case is the HCP crystal value. That is an ordered lattice, not a physical granular packing. A real monodisperse granular bed random-close-packs at about 0.64, and starting a simulation from a perfect crystal is usually undesirable anyway, since crystalline beds shear along lattice planes and give you the wrong dilatancy. So a polydisperse bed that settles to something below 0.74 is not failing. It is behaving correctly. With a sufficiently wide size distribution, small particles fill the interstices between large ones and you can comfortably exceed the monodisperse 0.64, which is the number worth comparing against.

 

Good luck with it,

Dan

------------------------------------------------

Robert and Laura Hensel Professor

NVIDIA CUDA Fellow

Department of Mechanical & Aerospace Engineering

Department of Electrical & Computer Engineering

Department of Computer Science

University of Wisconsin - Madison

4150ME, 1513 University Avenue

Madison, WI 53706-1572

608 772 0914

http://sbel.wisc.edu/

http://projectchrono.org/

------------------------------------------------

---

--
You received this message because you are subscribed to the Google Groups "ProjectChrono" group.
To unsubscribe from this group and stop receiving emails from it, send an email to projectchron...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/projectchrono/92e1b9ed-bbb3-44c3-8bfe-023879c83a4fn%40googlegroups.com.

[Shahriar Karim Shafin]

Thanks for the suggestion. 

I am encountering a problem after updating DEM-Engine. The same C++ simulation worked correctly with the previous version, but after updating to the current version, it fails with a “too many geometries in a bin” error during contact detection for all the codes with worked in previous version. I updated DEM-Engine because I had previously encountered the issue where the simulation would stop progressing while the GPU remained at 100% utilization, which was addressed by the recently merged fix I believe(#72; 1st GPU stuck at 100% and 2nd GPU 0% in the settling phase).

I reviewed the troubleshooting documentation and tried the recommended solutions, including performing clean rebuilds of both DEM-Engine and my application. However, the bin-related error persists. When I rebuild the application against the previous DEM-Engine version, the same simulation runs successfully again. The error -

-------- Simulation crashed "potentially" due to too many geometries in a bin --------
The dT reported max velocity is 0.000130473
------------------------------------
If the velocity is huge, then the simulation probably diverged due to encountering large particle velocities.
Decreasing the step size could help, and remember to check if your simulation objects are initially within the domain you specified.
------------------------------------
If the velocity is fair, and you *are* using a custom force model, one thing to do is to SetForceCalcThreadsPerBlock to a small number like 128 (see README.md troubleshooting for details).
If you are not using a custom model, one thing to do is to ensure the simulation world size (InstructBoxDomainDimension) is not orders of magnitude larger than the actual space the simulation entities take up.
If none works and you are going to discuss this on forum https://groups.google.com/g/projectchrono, please include a visual rendering of the simulation before crash.

terminate called after throwing an instance of 'std::runtime_error'
  what():  GPU Assertion: an illegal memory access was encountered. This happened in /data/lab/boddeti/DEM-Engine/src/algorithms/DEMCubContactDetection.cu:464  

My code (I have tried reducing Box dimesion too) -

   

          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.SetForceCalcThreadsPerBlock(128);

    //DEMSim.SetInitBinNumTarget(400000000);

    auto mat_type_container =

        DEMSim.LoadMaterial({{"E", 1e9}, {"nu", 0.3}, {"CoR", 0.1}, {"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);

    */

    // We can specify the force model using a file.

   

    auto my_force_model = DEMSim.ReadContactForceModel("ForceModelEEPA.cu");

    // Those following lines are needed. We must let the solver know that those var names are history variable etc.

    my_force_model->SetMustHaveMatProp({"E", "nu", "CoR", "mu", "Crr"});

    my_force_model->SetMustPairwiseMatProp({"CoR", "mu", "Crr"});

    // Pay attention to the extra per-contact wildcard `delta_max' here.

    my_force_model->SetPerContactWildcards(

        {"delta_time", "delta_tan_x", "delta_tan_y", "delta_tan_z", "delta_max"});

   

    path out_dir = current_path();

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

    create_directory(out_dir);  

   

    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 = 10.0;

    float container_diameter = 0.10;

    float terrain_density = 5.03e3;

   

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

    float sphere_rad = 0.0010;          // Mean Radius

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

    float compaction_speed_1 = -0.01f;

    float compaction_speed_2 = -0.002f;

    float compaction_speed_3 = 0.002f;  

    float step_size = 8e-7;

    //float fact_radius = 0.9;

    DEMSim.InstructBoxDomainDimension(world_size, world_size, world_size);

    DEMSim.InstructBoxDomainBoundingBC("none", mat_type_container);

    double bottom = 0;

    double top = 3.0;

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

    // ## 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.00, 0.00, 0 - 2 * sphere_rad);  // z

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

    walls->Scale(make_float3(4.0, 0.1, 12.0));

    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);

[Ruochun Zhang]

Hi Shahriar,


This is interesting... I would like to take a closer look. The code shown in the thread appears incomplete. Could you share with us the full code as attached files (including the custom force model ForceModelEEPA.cu)? Also, you mentioned that a previous version works, could you share that version's commit ID with me? Then I can start researching what went wrong.

Thank you,

Ruochun

[Shahriar Karim Shafin]

I have tried with the default Hertz-Mindlin Model also, same problem occurs, Currently I am using Edinburgh Elasto-Plastic Contact Model for plasticity and adhesion. The main problem with the previous version, sometimes the program freezes at 1 GPU utilization 100% (no progress) where the other gpu is idle at the start of settling. 

Previous Commit ID : 34d6f9b  (34d6f9baf3a290a2cf8526b6d81f06fb81b2f95d)

Commit message: : Fix a problem that could cause problem in submodule

 

[Ruochun Zhang]

Hi Shahriar,

It's actually running for me but the step size is small, so it might take a while before I see the problem. Could you tell me at which frame you get the error, generally? 

My experience is that if it fails with too many geos in bin, then either the bin size auto-adaptation failed, or the force model has no effect and you have gravity, so all particles overlap in the end, and a bin cannot host that many of them. So while I am doing that, that you could try is, first to use the newest main branch, then set a SetInitBinNumTarget (perhaps 1e7), then DisableAdaptiveBinSize(), to see if these can make your simulation run. If it still does not, we at least eliminate a possibility.

Ruochun

[Ruochun Zhang]

Hi Shahriar,

But I was able to run the simulation normally. I am using the default force model, and I attached a movie showing some 30 frames I got. Please have a look and see if it is what you'd expect.

Have you tried rebooting to clear the /tmp directory, and then doing a full clean rebuild? Normally this is not needed, but from what you are describing, I have a feeling the code is linking against some old JIT cache somehow. Not sure how this could happen. Trying it on another machine may also be helpful. And letting us know your GPU model can help too.

Since you said you can at least get frame 0 and 1, we should also have a look at those rendered frames, to make sure it represents normal physics. This is probably one of those system environment issues and might be a bit more painful...

Thank you,

Ruochun