Hi,
This is a DEME-related question.
I have been running into a problem where my simulation crashes after being
normal for a while. The error I get is the following:
//////////////////////////////////////////////////
-------- Simulation crashed potentially due to too many geometries in a bin
--------
Right now, the dT reported (by user specification or by calculation) max
velocity is 0.133465
The contact margin thickness is 9.35108e-06
If the velocity is extremely large, then the simulation probably diverged
due to encountering large particle velocities, and decreasing the step size
could help.
If the velocity is fair but the margin is large compared to particle sizes,
then perhaps too many contact geometries are in one bin, and decreasing the
step size, update frequency or the bin size could help.
If they are both fair and you do not see "exceeding maximum allowance"
reports before the crash, then it is probably not too many geometries in a
bin and it crashed for other reasons.
terminate called after throwing an instance of 'std::runtime_error'
what(): GPU Assertion: an illegal memory access was encountered. This
happened in /DEM-Engine/src/algorithms/DEMCubContactDetection.cu:
////////////////////////////////////////////////////////////////
I have tried to reduce my simulation bin size to as small as 0.5*particle
radius. I have also tried to reduce/increase other parameters, such as
update frequency and safety multiplier but still, the simulation crashes
after being normal for a while (I have a video that I could share with you
via email if you like). In addition, I have tried to reduce the time step
size very much (4e-7) but that did not seem to work. Also, I have reduced
my mesh to have a Total num of triangles: 6790 which I do not think is
really large. I have attached my sim file for your reference.
In addition, I have tried to use a different material from one of your
demos with the same time step but I still seem to have the same problem.
Also, one thing that I noticed, every time I increased the CDupdate
frequency value, the simulation reports a higher value of Average steps per
dynamic update. for example, when I set my update frequency to 15, the
simulation reports the Average steps per dynamic update: 16.94662. In
addition, when I increase my CD update to 20 the simulation reports the
Average steps per dynamic update: 21.997. Is that how it is supposed to be?
Thank you in advance for your help,
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/projectchrono/5bdfbb8e-9816-47d4-8f7a-15abf2f8a394n%40googlegroups.com.
// Copyright (c) 2021, SBEL GPU Development Team
// Copyright (c) 2021, University of Wisconsin - Madison
//
// SPDX-License-Identifier: BSD-3-Clause
#include <core/ApiVersion.h>
#include <core/utils/ThreadManager.h>
#include <DEM/API.h>
#include <DEM/HostSideHelpers.hpp>
#include <DEM/utils/Samplers.hpp>
#include <cstdio>
#include <chrono>
#include <filesystem>
using namespace deme;
using namespace std::filesystem;
int main() {
DEMSolver DEMSim;
DEMSim.SetVerbosity(INFO);
DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV);
DEMSim.SetOutputContent(OUTPUT_CONTENT::ABSV | OUTPUT_CONTENT:: VEL );
DEMSim.SetMeshOutputFormat(MESH_FORMAT::VTK);
DEMSim.EnsureKernelErrMsgLineNum();
// E, nu, CoR, mu, Crr...
auto mat_type_screw = DEMSim.LoadMaterial({{"E", 4.1e9}, {"nu", 0.394},
{"CoR", 0.3}, {"mu", 0.04}, {"Crr", 0.01}});
auto mat_type_terrain = DEMSim.LoadMaterial({{"E", 1e9}, {"nu", 0.45},
{"CoR", 0.5}, {"mu", 0.4}, {"Crr", 0.175}});
//auto mat_type_terrain = DEMSim.LoadMaterial({{"E", 1.0e9}, {"nu", 0.32},
{"CoR", 0.36}, {"mu", 0.5}, {"Crr", 1.00}});
// 1- Boundary:
float step_size = 5e-7;
//double world_size = 0.75;
double xBoundary = 0.35;
double yBoundary = 0.42;
double zBoundary = 0.38;
DEMSim.InstructBoxDomainDimension(xBoundary, yBoundary, zBoundary);
DEMSim.InstructBoxDomainBoundingBC("top_open", mat_type_terrain);
DEMSim.SetCoordSysOrigin("center");
// 2- Screw:
// At the beginning of the simulation and until the end of the copression,
the screw has a family type =1. this family is set to be fixed and does not
contact the particles.
// After, the screw family type is set to 3.
auto projectile = DEMSim.AddWavefrontMeshObject("./screwCM.obj",
mat_type_screw);
std::cout << "Total num of triangles: " << projectile->GetNumTriangles() <<
std::endl;
projectile->InformCentroidPrincipal(make_float3(0.095965,0.237612,0.095964),make_float4
( 0,0,0,1));
projectile->SetInitPos(make_float3(0, 0.1, 0.13 ));
float screw_mass = 0.59; // screw mass 0.59kg .
projectile->SetMass(screw_mass);
projectile->SetMOI(make_float3(0.062, 0.003, 0.0062));
projectile->SetFamily(1);
DEMSim.SetFamilyFixed(1);
DEMSim.DisableContactBetweenFamilies(0, 1);
// Track the projectile
auto proj_tracker = DEMSim.Track(projectile);
float rev_per_sec = 1;
float ang_vel_Z = rev_per_sec * 3.14;
float w_r = -1; //rad/s
// Screw Coredinate system with respect to global corrdinate system:
// y-aixs = z-axis global.
// z-aixs = y-axis global.
// x-aixs = x-axis global.
// In fact, because the cone's motion is completely pre-determined, we can
just prescribe family 3
//DEMSim.SetFamilyPrescribedLinVel(2,"4", "none", "none", false);
DEMSim.SetFamilyPrescribedLinVel(3, "0","none", "none", false);
DEMSim.AddFamilyPrescribedAcc(3, "none","none",
to_string_with_precision(-10.0/screw_mass));
DEMSim.SetFamilyPrescribedAngVel(3, "0", to_string_with_precision(w_r),
"0", false);
// 3-Particles:
float terrain_rad = 0.003;
auto template_terrain = DEMSim.LoadSphereType(terrain_rad * terrain_rad *
terrain_rad * 2.5e3 * 4.0 / 3.0 * 3.14,
terrain_rad,
mat_type_terrain);
std::vector<std::shared_ptr<DEMClumpTemplate>> clump_types;
clump_types.push_back(template_terrain);
// number if different clum types.
int num_template = 1;
// Generate initial clumps for piling
float spacing = 2.02*terrain_rad;
float sample_halfwidth_x = xBoundary / 2.0 - 2*terrain_rad;
float sample_halfwidth_y = yBoundary / 2.0 - 2*terrain_rad;
float fill_height = 0.2;//0.26;
float fill_bottom = -zBoundary / 2.0+ spacing;
PDSampler sampler(spacing);
// Use a PDSampler-based clump generation process
std::vector<std::shared_ptr<DEMClumpTemplate>> input_pile_template_type;
std::vector<float3> input_pile_xyz;
float layer_z = 0;
while (layer_z < fill_height) {
float3 sample_center = make_float3(0, 0, fill_bottom + layer_z );
auto layer_xyz = sampler.SampleBox(sample_center,
make_float3(sample_halfwidth_x, sample_halfwidth_y, 0));
unsigned int num_clumps = layer_xyz.size();
// Select from available clump types
for (unsigned int i = 0; i < num_clumps; i++) {
input_pile_template_type.push_back(clump_types.at(i %
num_template));
}
input_pile_xyz.insert(input_pile_xyz.end(), layer_xyz.begin(),
layer_xyz.end());
layer_z += spacing;
}
// Calling AddClumps a to add clumps to the system
auto the_pile = DEMSim.AddClumps(input_pile_template_type, input_pile_xyz);
std::cout << "Terrain loaded: " << std::endl;
size_t n_particles = input_pile_xyz.size();
std::cout << "Added Number of clump:" << n_particles << std::endl;
// Create a inspector to find out stuff
auto max_z_finder = DEMSim.CreateInspector("clump_max_z");
auto min_z_finder = DEMSim.CreateInspector("clump_min_z");
float max_z;
float min_z;
// 4- Compressor:
// Now add a plane to compress the sample
//auto compressor = DEMSim.AddExternalObject();
//compressor->AddPlane(make_float3(0, 0, 0), make_float3(0, 0, -1),
mat_type_terrain);
//compressor->SetFamily(10);
//DEMSim.DisableContactBetweenFamilies(0, 10);
//DEMSim.DisableContactBetweenFamilies(1, 10);
//DEMSim.DisableContactBetweenFamilies(2, 10);
//DEMSim.DisableContactBetweenFamilies(3, 10);
//DEMSim.SetFamilyFixed(10);
//auto compressor_tracker = DEMSim.Track(compressor);
// 5- Simulation initilization:
//DEMSim.SetCoordSysOrigin(make_float3(xBoundary/2,yBoundary/2,zBoundary/2));
DEMSim.SetInitTimeStep(step_size);
DEMSim.SetGravitationalAcceleration(make_float3(0, 0, -9.81));
// If you want to use a large UpdateFreq then you have to expand spheres to
ensure safety
DEMSim.SetCDUpdateFreq(15);
DEMSim.SetMaxVelocity(5.);
DEMSim.SetExpandSafetyMultiplier(1.1);
DEMSim.SetExpandSafetyAdder(1.1);
DEMSim.SetInitBinSize(2.*terrain_rad);
DEMSim.Initialize();
// 6- Drop Particles:
std::cout << "//////////////// Drop Particles //////////////////: " <<
std::endl;
path out_dir = current_path();
out_dir += "/ScrewDrop";
create_directory(out_dir);
float settling_end = 0.450;
float compression_end = 0.35;//1.050;
float rolling_end = 2.50;
unsigned int fps = 70;
float frame_time = 1.0 / fps;
float sim_time = 0.0;
unsigned int total_frames = (unsigned int)((rolling_end )* fps);
std::cout << "Output at " << fps << " FPS" << std::endl;
unsigned int currframe = 0;
for (float t = 0; t < settling_end; t += frame_time) {
std::cout << "Frame: " << currframe<< " of " << total_frames <<
std::endl;
char filename[200], meshfilename[200], cnt_filename[200];
sprintf(filename, "%s/DEMdemo_output_%04d.csv", out_dir.c_str(),
currframe);
sprintf(meshfilename, "%s/DEMdemo_mesh_%04d.vtk", out_dir.c_str(),
currframe);
DEMSim.WriteSphereFile(std::string(filename));
DEMSim.WriteMeshFile(std::string(meshfilename));
currframe++;
sim_time +=frame_time;
DEMSim.DoDynamics(frame_time);
DEMSim.ShowThreadCollaborationStats();
}
max_z = max_z_finder->GetValue();
min_z = min_z_finder->GetValue();
float porosity;
float total_volume_of_spheres =( 4.0 / 3.0) * 3.14159
*(terrain_rad*terrain_rad*terrain_rad)*n_particles;
float total_volume = (xBoundary * yBoundary ) * (max_z - min_z);
porosity = (total_volume-total_volume_of_spheres) / total_volume;
std::cout << "Max Z is: " << max_z << std::endl;
std::cout << "Min Z is: " << min_z << std::endl;
std::cout << "Porosity: " << porosity << std::endl;
// 7-Compresse:
//std::cout << "//////////////// Compression //////////////////: " <<
std::endl;
//DEMSim.EnableContactBetweenFamilies(0, 10);
//double compressor_vel = 0.1;
//float terrain_max_z = max_z_finder->GetValue();
//compressor_tracker->SetPos(make_float3(0, 0, terrain_max_z));
// for (float t = sim_time; t < compression_end; t += frame_time) {
// std::cout << "Frame: " << currframe<< " of " << total_frames <<
std::endl;
// char filename[200], meshfilename[200], cnt_filename[200];
// sprintf(filename, "%s/DEMdemo_output_%04d.csv", out_dir.c_str(),
currframe);
// sprintf(meshfilename, "%s/DEMdemo_mesh_%04d.vtk", out_dir.c_str(),
currframe);
// DEMSim.WriteSphereFile(std::string(filename));
// DEMSim.WriteMeshFile(std::string(meshfilename));
// currframe++;
//terrain_max_z -= compressor_vel * frame_time;
//compressor_tracker->SetPos(make_float3(0, 0, terrain_max_z));
// proj_tracker->SetPos(make_float3(0, 0, terrain_max_z+0.05)); //0.2
screw diamerter
// DEMSim.DoDynamicsThenSync(frame_time);
// DEMSim.ShowThreadCollaborationStats();
// total_volume_of_spheres =( 4.0 / 3.0) * 3.14159
*(terrain_rad*terrain_rad*terrain_rad)*n_particles;
// total_volume = (xBoundary * yBoundary ) * (max_z - min_z);
// porosity = (total_volume-total_volume_of_spheres) / total_volume;
// max_z = max_z_finder->GetValue();
// min_z = min_z_finder->GetValue();
// std::cout << "Max Z is: " << max_z << std::endl;
// std::cout << "Min Z is: " << min_z << std::endl;
// std::cout << "Porosity: " << porosity << std::endl;
// sim_time +=frame_time;
// std::cout << "simulation time: " << sim_time << std::endl;
//}
//Remove compressor
//DEMSim.DisableContactBetweenFamilies(0, 10);
//DEMSim.DisableContactBetweenFamilies(1, 10);
//DEMSim.DisableContactBetweenFamilies(2, 10);
//DEMSim.DisableContactBetweenFamilies(3, 10);
DEMSim.DoDynamicsThenSync(frame_time);
std::cout << " Screw Postion Set " << std::endl;
float terrain_max_z = max_z_finder->GetValue();
proj_tracker->SetPos(make_float3(0, 0.1, terrain_max_z+0.1)); //0.2 screw
diamerter
DEMSim.DoDynamicsThenSync(frame_time);
DEMSim.ShowThreadCollaborationStats();
// Drop Screw:
DEMSim.ChangeFamily(1, 3);
std::cout << "//////////////// Spen Screw //////////////////: " <<
std::endl;
for (float t = sim_time; t < rolling_end; t += frame_time) {
std::cout << "Frame: " << currframe<< " of " << total_frames <<
std::endl;
char filename[200], meshfilename[200], cnt_filename[200];
sprintf(filename, "%s/DEMdemo_output_%04d.csv", out_dir.c_str(),
currframe);
sprintf(meshfilename, "%s/DEMdemo_mesh_%04d.vtk", out_dir.c_str(),
currframe);
DEMSim.WriteSphereFile(std::string(filename));
DEMSim.WriteMeshFile(std::string(meshfilename));
float3 pos_screw = proj_tracker->Pos();
float3 force = proj_tracker->ContactAcc();
std::cout << "screw postion: " << pos_screw.x << ", " << pos_screw.y <<
", " << pos_screw.z << std::endl;
std::cout << "screw forces: " << force.x << ", " << force.y << ", " <<
force.z << std::endl;
currframe++;
sim_time +=frame_time;
std::cout << " --Sync-- " << std::endl;
DEMSim.DoDynamics(frame_time);
DEMSim.ShowThreadCollaborationStats();
std::cout << " --Sync Complete-- " << std::endl;
}
std::cout << "simulation time: " << sim_time << std::endl;
std::cout << "DEMdemo_ScrewDrop exiting..." << std::endl;
return 0;
}
