Particle spawning issues

[Liam Murray]

I have tried to make progress on this simulation from the last time i posted on here. I got the shell to spin successfully and the meshes lined up thanks to feedback. What im trying to do is fill up the basin and shell meshes sort of like the mixer demo so that it is filled all of the way when the simulation initializes to start getting data and not having to wait for the particles to convey to the top. Ive tried messing with a lot of these parameters, but im just continuing to get errors like " outputs this error: I know setting max sphere in bin isnt the greatest practice or max contacts, but i keep running into this again and again and i wish i could spawn it in cleanly into the basin. the image attached is what the setup looks like without particles currently.


 butterco@NASALaptop:~/DEM-Engine/build/bin$ ./augertest Total number of particles generated: 71996 WARNING! One GPU device is detected. On consumer cards, DEME's performance edge is limited with only one GPU. Try allocating 2 GPU devices if possible. WARNING! Some clumps do not have their family numbers specified, so defaulted to 0 Number of total active devices: 1 User-specified X-dimension range: [-0.75, 0.75] User-specified Y-dimension range: [-0.75, 0.75] User-specified Z-dimension range: [-1, 1] User-specified dimensions should NOT be larger than the following simulation world. The dimension of the simulation world: 1.7999999523162842, 1.7999999523162842, 3.5999999046325684 Simulation world X range: [-0.9, 0.9] Simulation world Y range: [-0.9, 0.9] Simulation world Z range: [-1.2, 2.4] The length unit in this simulation is: 1.3096723358585471e-11 The edge length of a voxel: 8.5830686202825746e-07 The initial time step size: 5e-06 The initial edge length of a bin: 0.021110625075200021 The initial number of bins: 1264716 The total number of clumps: 71996 The combined number of component spheres: 215988 The total number of analytical objects: 1 The total number of meshes: 3 Grand total number of owners: 72000 The number of material types: 1 History-based Hertzian contact model is in use. The solver to set to adaptively change the contact margin size. Bin 884 contains 35373 sphere components, exceeding maximum allowance (32768). If you want the solver to run despite this, set allowance higher via SetMaxSphereInBin before simulation starts. -------- Simulation crashed "potentially" due to too many geometries in a bin -------- The dT reported max velocity is 0 ------------------------------------ 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 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: unspecified launch failure. This happened in /home/butterco/DEM-Engine/src/algorithms/DEMCubContactDetection.cu:384 Aborted (core dumped) butterco@NASALaptop:~/DEM-Engine/build/bin$"

"

#include <DEM/API.h>
#include <DEM/HostSideHelpers.hpp>
#include <DEM/utils/Samplers.hpp>
#include <chrono>
#include <iostream>
#include <filesystem>

using namespace deme;
using namespace std::chrono;
namespace fs = std::filesystem;

int main() {
    DEMSolver DEMSim;
    DEMSim.SetVerbosity(INFO);
    DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV);
    DEMSim.SetOutputContent(OUTPUT_CONTENT::ABSV);
    DEMSim.SetMeshOutputFormat(MESH_FORMAT::VTK);

    // Material properties
    auto mat_type_mixer = DEMSim.LoadMaterial({{"E", 5e5}, {"nu", 0.3}, {"CoR", 0.2}, {"mu", 0.3}, {"Crr", 0.0}});

    // Set simulation domain dimensions
    DEMSim.InstructBoxDomainDimension(1.5, 1.5, 2.0);
    DEMSim.InstructBoxDomainBoundingBC("all", mat_type_mixer);

    // Load mesh objects (Basin, Auger, Drum)
    auto basin = DEMSim.AddWavefrontMeshObject("/home/butterco/DEM-Engine/data/mesh/augertest_shapes/body_1_1_collision.obj", mat_type_mixer);
    auto auger = DEMSim.AddWavefrontMeshObject("/home/butterco/DEM-Engine/data/mesh/augertest_shapes/body_2_1_collision.obj", mat_type_mixer);
    auto drum = DEMSim.AddWavefrontMeshObject("/home/butterco/DEM-Engine/data/mesh/augertest_shapes/body_3_1_collision.obj", mat_type_mixer);

    // Set initial positions and orientations
    basin->Move(make_float3(-0.1, -0.1, -0.70), make_float4(0, 0, 0, 1));  
    drum->Move(make_float3(0, 0, -0.6), make_float4(0, 0.7071, 0.7071, 0));  
    auger->Move(make_float3(0, 0, 0.85), make_float4(1, 0, 0, 0));  

    // Set family numbers for simulation properties
    basin->SetFamily(0);  
    auger->SetFamily(1);  
    drum->SetFamily(2);  
    DEMSim.SetFamilyFixed(0);  
    DEMSim.SetFamilyFixed(1);  
    DEMSim.SetFamilyPrescribedAngVel(2, "0", "0", "0.1");

    // Define particle template with increased radius
    float particle_radius = 0.005f; // Increased radius to reduce particle count
    float particle_volume = (4.0f / 3.0f) * 3.14159f * std::pow(particle_radius, 3);
    float particle_mass = 2.6e3f * particle_volume;  // Density * volume (kg)
    float3 particle_moi = make_float3(1.8327927f, 2.1580013f, 0.77010059f) * (std::pow(particle_radius, 2) / std::pow(0.005f, 2));

    auto particle_template = DEMSim.LoadClumpType(particle_mass, particle_moi, "/home/butterco/DEM-Engine/data/clumps/3_clump.csv", mat_type_mixer);
    // No scaling needed as the particle radius matches the clump file

    // Sampler setup with increased spacing
    HCPSampler basin_sampler(3.0f * particle_radius); // Increased spacing
    float3 basin_center = make_float3(0.0f, 0.0f, -0.65f);
    float basin_radius = 0.07f;
    float basin_height = 0.15f;
    auto basin_positions = basin_sampler.SampleCylinderZ(basin_center, basin_radius, basin_height);

    std::cout << "Total number of particles generated: " << basin_positions.size() << std::endl;

    if (!basin_positions.empty()) {
        DEMSim.AddClumps(particle_template, basin_positions);
    } else {
        std::cerr << "Error: No particles were generated. Adjust sampling parameters." << std::endl;
        return -1;
    }

    DEMSim.SetInitTimeStep(5e-6f);
    DEMSim.SetGravitationalAcceleration(make_float3(0.0f, 0.0f, -9.81f));
    DEMSim.SetCDUpdateFreq(40);
    DEMSim.SetExpandSafetyAdder(2.0f);
    DEMSim.SetErrorOutAvgContacts(50);
    DEMSim.SetAdaptiveBinSizeUpperProactivity(0.5f);
    DEMSim.SetAdaptiveBinSizeLowerProactivity(0.15f);

    // Set initial bin size to prevent bins from becoming overpopulated
    DEMSim.SetInitBinSize(25.0f * particle_radius);

    // Increase maximum spheres per bin if necessary
    DEMSim.SetMaxSphereInBin(50000);

    // Initialize the simulation
    DEMSim.Initialize();

    // Create output directory
    std::string output_dir = "/home/butterco/DEM-Engine/data/output/";
    fs::create_directories(output_dir);

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

    // Set simulation time and frame rate
    float sim_end = 1.0f;  
    unsigned int fps = 10;  
    float frame_time = 1.0f / fps;
    unsigned int curr_frame = 0;

    // Simulation loop
    for (float t = 0.0f; t < sim_end; t += frame_time) {
        std::cout << "Simulating frame " << curr_frame << std::endl;

        char filename[200], mesh_filename[200];
        sprintf(filename, "%sDEM_output_%04d.csv", output_dir.c_str(), curr_frame);
        sprintf(mesh_filename, "%sDEM_mesh_output_%04d.vtk", output_dir.c_str(), curr_frame++);
        DEMSim.WriteSphereFile(std::string(filename));
        DEMSim.WriteMeshFile(std::string(mesh_filename));

        DEMSim.DoDynamics(frame_time);

        float max_velocity = max_velocity_finder->GetValue();
        std::cout << "Max velocity of any particle: " << max_velocity << std::endl;
    }

    DEMSim.ShowTimingStats();
    std::cout << "Simulation complete." << std::endl;

    return 0;
}

"

[Liam Murray]

the thing im having trouble seeing possible is the mixer demo spawns the particles above the mesh when the demo initializes and im not sure compared to programs like StarCCM+ or ansys Rocky where there are options for assisting in it not spawning within the mesh. would it be something of the sort of making a solver within the code given the collision file to have the particles not contact the mesh, but im not sure how one would code that exactly, im just putting my thoughts out there. or is it better to spawn a lot of particles spawn in above the mesh and wait for all of them to settle  

[Ruochun Zhang]

Hi Liam,

We see your simulation crashed at the first step, which shows problems with the initial setup. 

In your script, you didn't scale the clump info in the 3_clump.csv file, and left a note saying there's no need. Did you modify the file or is it still the same file bundled with DEM-Engine? If it's the same, then the clump has a size ≈ 2, as large as your simulation domain. This can explain the error you had, since if 70000 such particles overlap together, no matter how you divide the domain, a bin will intersect with most of them.

In the picture you did not render the particles, but why? You need to visualize them to understand the problem.

About the utility that prevents particles from spawning in the mesh, I think it's a nice-to-have and won't be too difficult to implement, but it is currently not our priority to make it a part of the package. You could easily write it yourself using C++ or Python though. DEM-Engine has support for smartly disabling these particles before the simulation starts, but it's an advanced usage. For now I suggest sampling the particles at safe locations and then letting them settle.

Thank you,

Ruochun 

[Liam Murray]

Thanks for the response! 

I updated the code to this as reccomended and ran a lot of different scenarios and began to think that maybe I don't need the basin and could have a particle domain that just surrounds the shell as the image attached shows. I've been really trying to show some decent progress since we've last talked. Now I just need to densely pack the particles and decrease the size eventually to start conveying some stuff! Then maybe trying to write the code for preventing it spawning in mesh. Let me know what you think

#include <DEM/API.h>

#include <DEM/HostSideHelpers.hpp>

#include <DEM/utils/Samplers.hpp>

#include <chrono>

#include <iostream>

#include <filesystem>

#include <cmath> // For sqrt function

using namespace deme;

using namespace std::chrono;

namespace fs = std::filesystem;

int main() {

    // Initialize the DEM solver

    DEMSolver DEMSim;

    DEMSim.SetVerbosity(INFO);

    DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV);

    DEMSim.SetOutputContent(OUTPUT_CONTENT::ABSV);

    DEMSim.SetMeshOutputFormat(MESH_FORMAT::VTK);

    // Material properties with reduced stiffness and increased damping

    auto mat_type_mixer = DEMSim.LoadMaterial({

        {"E", 1e5f},    // Reduced stiffness to prevent high contact forces

        {"nu", 0.3f},

        {"CoR", 0.2f},

        {"mu", 0.3f},

        {"Crr", 0.2f}    // Increased rolling resistance

    });

    // Adjusted simulation domain dimensions

    float domain_x = 3.0f; // Increased domain to accommodate spawning area

    float domain_y = 3.0f;

    float domain_z = 3.0f;

    DEMSim.InstructBoxDomainDimension(domain_x, domain_y, domain_z);

    DEMSim.InstructBoxDomainBoundingBC("none", mat_type_mixer); // We'll add custom walls instead

    // Load mesh objects (Auger and Drum) without scaling

    // Commented out the basin as per user request

    /*

    auto basin = DEMSim.AddWavefrontMeshObject(

        "/home/butterco/DEM-Engine/data/mesh/augertest_shapes/body_1_1_collision.obj",

        mat_type_mixer

    );

    basin->SetFamily(0);

    DEMSim.SetFamilyFixed(0); // Basin is fixed

    */

    auto auger = DEMSim.AddWavefrontMeshObject(

        "/home/butterco/DEM-Engine/data/mesh/augertest_shapes/body_2_1_collision.obj",

        mat_type_mixer

    );

    auto drum = DEMSim.AddWavefrontMeshObject(

        "/home/butterco/DEM-Engine/data/mesh/augertest_shapes/body_3_1_collision.obj",

        mat_type_mixer

    );

   

    drum->Move(make_float3(0, 0, -0.6f), make_float4(0, 0.7071f, 0.7071f, 0));

    auger->Move(make_float3(0, 0, 0.85f), make_float4(1, 0, 0, 0));

    // Set family numbers for simulation properties

    // basin->SetFamily(0); // Basin is commented out

    auger->SetFamily(1);

    drum->SetFamily(2);

    // Fix the auger

    DEMSim.SetFamilyFixed(1); // Auger is completely fixed

   

    // DEMSim.SetFamilyFixed(2); // Removed to allow drum rotation

    // Assign prescribed angular velocity to the drum

    DEMSim.SetFamilyPrescribedAngVel(2, "0", "0", "0.5");  // Drum rotates around Z-axis at 0.5 rad/s

    // Define particle template with decreased radius

    float particle_radius = 0.01f; // Particle radius decreased to 10 mm

    float mass = 2.6e3f * (4.0f / 3.0f) * M_PI * pow(particle_radius, 3); // Mass in kg (~0.1088 kg)

    float3 MOI = (2.0f / 5.0f) * mass * particle_radius * particle_radius * make_float3(1.0f, 1.0f, 1.0f); // MOI for a sphere (~4.352e-5 kg·m²)

    auto particle_template = DEMSim.LoadSphereType(particle_radius, mass, mat_type_mixer);

    // Sampler setup to spawn particles around the drum and auger

    float drum_radius = 0.02f; // Drum radius in meters

    // Reverted inner and outer radii

    float inner_radius = 0.036f; // Original inner radius (3.6 cm)

    float outer_radius = 0.15f;  // Reduced outer radius to 0.15 meters to match particle domain

    // Adjusted fill height and fill bottom based on spawning volume calculation

    float fill_height = 0.2f; // Fill height in meters (20 cm)

    float fill_bottom = -0.7f; // Bottom of the fill region in meters

    float3 fill_center = make_float3(0.0f, 0.0f, fill_bottom + fill_height / 2.0f); // (0, 0, -0.6)

    // Debugging: Print inner and outer radii and fill height

    std::cout << "Inner Radius: " << inner_radius << " meters" << std::endl;

    std::cout << "Outer Radius: " << outer_radius << " meters" << std::endl;

    std::cout << "Fill Height: " << fill_height << " meters" << std::endl;

    // Adjust sampler spacing to prevent overlaps and slightly increase particle count

    HCPSampler sampler(2.8f * particle_radius); // Sampler spacing set to 0.028m (2.8x particle_radius)

    // Sample particles in the cylinder with outer radius

    auto input_xyz_all = sampler.SampleCylinderZ(fill_center, outer_radius, fill_height);

    // Filter out particles within the inner radius to create a donut shape

    std::vector<float3> input_xyz;

    for (const auto& pos : input_xyz_all) {

        float dx = pos.x - fill_center.x;

        float dy = pos.y - fill_center.y;

        float dist = sqrt(dx * dx + dy * dy);

        if (dist >= inner_radius && dist <= outer_radius) {

            input_xyz.push_back(pos);

        }

    }

    std::cout << "Total number of particles generated: " << input_xyz.size() << std::endl;

    if (!input_xyz.empty()) {

        // Add clumps and get the clump batch

        auto clump_batch = DEMSim.AddClumps(particle_template, input_xyz);

        clump_batch->SetFamily(3); // Assign particles to family 3, which is not fixed

    } else {

        std::cerr << "Error: No particles were generated. Adjust sampling parameters." << std::endl;

        return -1;

    }

    // Adjust the initial time step size

    DEMSim.SetInitTimeStep(0.5e-6f); // Decrease time step for higher stability

    // Set gravity and other simulation parameters

    DEMSim.SetGravitationalAcceleration(make_float3(0.0f, 0.0f, -9.81f));

    DEMSim.SetCDUpdateFreq(40);

    DEMSim.SetExpandSafetyAdder(2.0f);

    // Increase the maximum allowed average contacts per sphere

    DEMSim.SetErrorOutAvgContacts(300); // Increased from 250 to 300

    // Increase the allowed system max velocity

    DEMSim.SetErrorOutVelocity(3000.0f); // Increased from 2000 to 3000

    DEMSim.SetAdaptiveBinSizeUpperProactivity(0.5f);

    DEMSim.SetAdaptiveBinSizeLowerProactivity(0.15f);

    // Set initial bin size to prevent bins from becoming overpopulated

    DEMSim.SetInitBinSize(25.0f * particle_radius); // 25 * 0.01 = 0.25m

    // Increase maximum spheres per bin if necessary

    DEMSim.SetMaxSphereInBin(50000);

    // Add walls to prevent particles from falling over the edge

    auto walls = DEMSim.AddExternalObject();

    walls->AddPlane(make_float3(0, 0, -0.7f), make_float3(0, 0, 1), mat_type_mixer); // Bottom plane at Z=-0.7m

    walls->AddCylinder(

        make_float3(0, 0, 0),

        make_float3(0, 0, 1),

        outer_radius, // 0.15m

        mat_type_mixer

    ); // Outer cylinder

    walls->SetFamily(4); // Assign walls to a new family

    DEMSim.SetFamilyFixed(4); // Fix the walls

    // Initialize the simulation

    DEMSim.Initialize();

    // Create output directory

    std::string output_dir = "/home/butterco/DEM-Engine/data/output/";

    fs::create_directories(output_dir);

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

    // Extend the simulation duration

    float sim_end = 20.0f; // Simulation end time in seconds

    unsigned int fps = 20; // Frames per second

    float frame_time = 1.0f / fps;

    unsigned int curr_frame = 0;

    // Simulation loop

    for (float t = 0.0f; t < sim_end; t += frame_time) {

        std::cout << "Simulating frame " << curr_frame << std::endl;

        char filename[200], mesh_filename[200];

        sprintf(filename, "%sDEM_output_%04d.csv", output_dir.c_str(), curr_frame);

        sprintf(mesh_filename, "%sDEM_mesh_output_%04d.vtk", output_dir.c_str(), curr_frame++);

        DEMSim.WriteSphereFile(std::string(filename));

        DEMSim.WriteMeshFile(std::string(mesh_filename));

        DEMSim.DoDynamics(frame_time);

        float max_velocity = max_velocity_finder->GetValue();

        std::cout << "Max velocity of any particle: " << max_velocity << std::endl;

    }

    DEMSim.ShowTimingStats();

    std::cout << "Simulation complete." << std::endl;

    return 0;

}

best,

Liam

--
You received this message because you are subscribed to a topic in the Google Groups "ProjectChrono" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/projectchrono/ewhC4AYLF-Y/unsubscribe.
To unsubscribe from this group and all its topics, send an email to projectchron...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/projectchrono/c7552937-946a-4488-a32c-27cc429f24c9n%40googlegroups.com.

[Ruochun Zhang]

Hi Liam,

You know your needs the best, how you do it depends on you.

But again, to get a densely packed configuration, you'll usually likely need to run the simulation and let it settle. Looks like from the picture you still have initial mesh-particle penetration, so running it probably won't go well. For me, it won't be difficult to carefully select a couple of regions in the space so after spawning the particles, there won't be initial penetration; or, I feel the easiest way is still sampling a big box region without the mesh, then add the mesh into the simulation scene, even if it involves simulating sticking the mesh into the pile of particles, which is not too difficult.

Thank you,

Ruochun

[Liam Murray]

We got it going! Thank you again for all your continued help. Wanted to show you some good progress 

To view this discussion visit https://groups.google.com/d/msgid/projectchrono/33ae610f-5c39-4263-8f59-24b464cae99en%40googlegroups.com.

[Liam Murray]

best,

Liam Murray
AeroFly LLC
D1 Swimmer
SDSU

[Ruochun Zhang]

That's good to know. From your rendering, particles appear to be floating and sphere components appear disconnected, which indicates that the rendered radius is smaller than the actual number. Remember in Paraview, you need to scale the "r" column by 2 in the Glyph setting.

Thank you,

Ruochun

[Liam Murray]

It does look a heck of a lot better I had r=1 thank you! I do appreciate again all the time you've taken to help me get this far

best,

Liam Murray
D1 Swimmer
SDSU

To view this discussion visit https://groups.google.com/d/msgid/projectchrono/3405bac5-e865-41a2-8516-60a76a30a782n%40googlegroups.com.