Hi,
This is a DEME-related question.
I have been running a simulation for settling a cylindrical bed of
particles with a specific set of parameters. I have successfully been able
to settle the bed. However, when I try to change the "Young modulus" for
the particles from 1e9Pa to 1e8Pa, the simulation gives me an error that
says the bin contains ### of particles that is larger than 256 maxima. I
have tried to reduce my ben size using *SetInitBinSize. *However, the
simulation gives me the same error until I set the bin size to *1.55*particle
radius* at which the error message changes to "The simulation world has
4301703661 bins (for domain partitioning in contact detection), but the
largest bin ID that we can have is 4294967295." Please note that for a
Young modulus of 1e9Pa, the bin size was set to be a *2*particle radius*
and the simulation worked with no problem. The simulation file is attached
for your reference
Thank you so much in advance,
--
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/2db59f80-2864-4d0f-ade0-d81a16888a72n%40googlegroups.com.
//Cone Penetraion test: Mohammad
#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);
DEMSim.SetMeshOutputFormat(MESH_FORMAT::VTK);
DEMSim.SetContactOutputContent(OWNER | FORCE | POINT);
// Material Properties for cone and particles:
// E, nu, CoR, mu, Crr, CED...
auto mat_type_cone = DEMSim.LoadMaterial({{"E", 2.0e11}, {"nu", 0.27},
{"CoR", 0.88}, {"mu", 0.5}, {"Crr", 1.00}});
auto mat_type_particle = DEMSim.LoadMaterial({{"E", 1.0e8}, {"nu", 0.32},
{"CoR", 0.36}, {"mu", 0.5}, {"Crr", 1.00},{"CED", 20.0e6}});
auto mat_type_container = DEMSim.LoadMaterial({{"E", 1.0e8}, {"nu",
0.32}, {"CoR", 0.36}, {"mu", 0.5}, {"Crr", 1.00}});
// Cone & cylindrical bin Parameters:
double time_settle =0.5;
double time_compress = 0.65;
double time_decompress = 0.75;
//float cone_speed = 0.1;//m/s
float step_size = 5.0e-7;//s
double world_size = 4.0;
double soil_bin_diameter = 0.1; //m
//double cone_diameter = 0.01;//m
//float shaft_mass = 47.436 ;
//float tip_mass = 1.8137 ;
//double starting_height = -0.02;
float particle_radius = 0.002; //m
float sphere_density = 920.0;//kg/m3
// 1- Particles Initiation:
auto template_sphere = DEMSim.LoadSphereType(particle_radius *
particle_radius * particle_radius * sphere_density * 4 / 3 * 3.14159,
particle_radius,
mat_type_particle);
std::vector<std::shared_ptr<DEMClumpTemplate>> clump_types;
clump_types.push_back(template_sphere);
// Generate initial clumps for piling
int num_template = 1;
// Generate initial clumps for piling
float spacing = 2.00*particle_radius;
float fill_width = soil_bin_diameter / 2.0 - particle_radius;
float fill_height = 0.44;
float fill_bottom = -0.01 + spacing;
PDSampler sampler(spacing);
//HCPSampler 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);//+
spacing / 10
auto layer_xyz = sampler.SampleCylinderZ(sample_center, fill_width, 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);
size_t n_particles = input_pile_xyz.size();
std::cout << "Added clump type " << 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;
// 2- Domain Initiation:
DEMSim.InstructBoxDomainDimension(world_size, world_size, world_size);
// No need to add simulation `world' boundaries, b/c we'll add a
cylinderical container manually
DEMSim.InstructBoxDomainBoundingBC("none", mat_type_container);
// Now add a cylinderical boundary along with a bottom plane
auto walls = DEMSim.AddExternalObject();
walls->AddCylinder(make_float3(0), make_float3(0, 0, 1), soil_bin_diameter
/ 2.0, mat_type_container, 0);
walls->AddPlane(make_float3(0, 0, -0.01), make_float3(0, 0, 1),
mat_type_container);
// 3- 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_container);
compressor->SetFamily(10);
DEMSim.DisableContactBetweenFamilies(0, 10);
DEMSim.SetFamilyFixed(10);
auto compressor_tracker = DEMSim.Track(compressor);
// A custom force model can be read in through a file and used by the
simulation. Magic, right?
auto my_force_model = DEMSim.ReadContactForceModel("SJKR-E.cu");
// This custom force model still uses contact history arrays, so let's
define it
//my_force_model->SetPerContactWildcards({"delta_tan_x", "delta_tan_y",
"delta_tan_z"});
my_force_model->SetPerContactWildcards({"delta_tan_x", "delta_tan_y",
"delta_tan_z", "delta_time"});
DEMSim.SetCoordSysOrigin("center");
// 3- Simulation Parameters & initilization:
DEMSim.SetInitTimeStep(step_size);
DEMSim.SetGravitationalAcceleration(make_float3(0, 0, -9.81));
DEMSim.SetCDUpdateFreq(6);
DEMSim.SetMaxVelocity(25.);
DEMSim.SetExpandSafetyMultiplier(2.1);
DEMSim.SetInitBinSize(1.55* particle_radius);
DEMSim.SetIntegrator(TIME_INTEGRATOR::CENTERED_DIFFERENCE);
DEMSim.Initialize();
DEMSim.UpdateSimParams();
// 4-Runing Simulation:
path out_dir = current_path();
out_dir += "/DemoOutput_Pile";
create_directory(out_dir);
float curr_time_of_sim = 0; // keeps track of simulation time
unsigned int currframe = 0;
float out_fps = 100.0;
unsigned int out_steps = (unsigned int)(1 / (out_fps * step_size));
unsigned int total_frames = (unsigned int)(time_decompress * out_fps);
unsigned int curr_step = 0;
for (double t = 0; t < (double)time_settle; t += step_size, curr_step++) {
if (curr_step % out_steps == 0) {
char filename[100];
sprintf(filename, "%s/step%06d.csv", out_dir.c_str(), currframe++);
DEMSim.WriteSphereFile(std::string(filename));
std::cout << "Output frame " << currframe + 1 << " of " << total_frames
<< std::endl;
}
DEMSim.DoDynamics(step_size);
curr_time_of_sim += step_size;
}
auto total_volume_finder = DEMSim.CreateInspector("clump_volume");
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
*(particle_radius*particle_radius*particle_radius)*n_particles;
float calc =( 4 / 3) * 3.14159
*(particle_radius*particle_radius*particle_radius);
float total_volume = 3.14159 * (soil_bin_diameter * soil_bin_diameter /
4.0) * (max_z - min_z);
porosity = (total_volume-total_volume_of_spheres) / total_volume;
std::cout << "total_volume_of_spheres is: " << total_volume_of_spheres <<
std::endl;
std::cout << "calc total_volume_of_spheres is: " << calc << std::endl;
std::cout << "n_particles is: " << n_particles << std::endl;
std::cout << "total_volume is: " << total_volume << std::endl;
std::cout << "Max Z is: " << max_z << std::endl;
std::cout << "Min Z is: " << min_z << std::endl;
std::cout << "Porosity Before Compression: " << porosity << std::endl;
std::cout << "curr_time_of_sim1: " << curr_time_of_sim << std::endl;
DEMSim.EnableContactBetweenFamilies(0, 10);
//Compresse
std::cout << "//////////////// Compression //////////////////: " <<
std::endl;
std::cout << "here1: " << std::endl;
double compressor_vel = 0.1;
float terrain_max_z = max_z_finder->GetValue();
float c_max_z = max_z_finder->GetValue();
compressor_tracker->SetPos(make_float3(0, 0, terrain_max_z));
curr_step = 0;
for (double t = 0; t < (double)time_compress; t += step_size, curr_step++) {
float matter_volume = ( 4.0 / 3.0) * 3.14159
*(particle_radius*particle_radius*particle_radius)*n_particles;
float total_volume2 = 3.14159 * (soil_bin_diameter * soil_bin_diameter
/ 4.0) * (c_max_z - min_z);
double porosity = (total_volume2-matter_volume) / total_volume2;
float3 Position_of_compressor = compressor_tracker->Pos();
if (curr_step % out_steps == 0) {
char filename[100];
sprintf(filename, "%s/step%06d.csv", out_dir.c_str(), currframe++);
DEMSim.WriteSphereFile(std::string(filename));
std::cout << "Output frame " << currframe + 1 << " of " <<
total_frames << std::endl;
std::cout << "Compression porosity: " << porosity << std::endl;
std::cout << "Compression particle_max_z: " << c_max_z << std::endl;
}
terrain_max_z -= compressor_vel * step_size;
compressor_tracker->SetPos(make_float3(0, 0, terrain_max_z));
c_max_z = max_z_finder->GetValue();
DEMSim.DoDynamics(step_size);
curr_time_of_sim += step_size;
}
std::cout << "curr_time_of_sim2: " << curr_time_of_sim << std::endl;
std::cout << "//////////////// DeCompression //////////////////: " <<
std::endl;
double decompressor_vel = 0.01;
c_max_z = max_z_finder->GetValue();
terrain_max_z = max_z_finder->GetValue();
compressor_tracker->SetPos(make_float3(0, 0, terrain_max_z));
for (double t = 0; t < (double)time_decompress; t += step_size,
curr_step++) {
float matter_volume = ( 4.0 / 3.0) * 3.14159
*(particle_radius*particle_radius*particle_radius)*n_particles;
float total_volume2 = 3.14159 * (soil_bin_diameter * soil_bin_diameter
/ 4.0) * (c_max_z - min_z);
double porosity = (total_volume2-matter_volume) / total_volume2;
float3 Position_of_compressor = compressor_tracker->Pos();
if (curr_step % out_steps == 0) {
char filename[100];
sprintf(filename, "%s/step%06d.csv", out_dir.c_str(), currframe++);
DEMSim.WriteSphereFile(std::string(filename));
std::cout << "Output frame " << currframe + 1 << " of " <<
total_frames << std::endl;
std::cout << "Decompression porosity: " << porosity << std::endl;
std::cout << "Decompression particle_max_z: " << c_max_z <<
std::endl;
}
terrain_max_z += decompressor_vel * step_size;
compressor_tracker->SetPos(make_float3(0, 0, terrain_max_z));
c_max_z = max_z_finder->GetValue();
DEMSim.DoDynamics(step_size);
curr_time_of_sim += step_size;
}
std::cout << "curr_time_of_sim3: " << curr_time_of_sim << std::endl;
//Remove compressor
DEMSim.DoDynamicsThenSync(step_size);
DEMSim.DisableContactBetweenFamilies(0, 10);
DEMSim.DoDynamicsThenSync(step_size);
terrain_max_z = max_z_finder->GetValue();
for (double t = 0; t < 0.05; t += step_size, curr_step++) {
float matter_volume = ( 4.0/ 3.0) * 3.14159
*(particle_radius*particle_radius*particle_radius)*n_particles;
float total_volume2 = 3.14159 * (soil_bin_diameter * soil_bin_diameter
/ 4.0) * (c_max_z - min_z);
double porosity = (total_volume2-matter_volume) / total_volume2;
float3 Position_of_compressor = compressor_tracker->Pos();
if (curr_step % out_steps == 0) {
char filename[100];
sprintf(filename, "%s/step%06d.csv", out_dir.c_str(), currframe++);
DEMSim.WriteSphereFile(std::string(filename));
std::cout << "Output frame " << currframe + 1 << " of " <<
total_frames << std::endl;
std::cout << " final porosity: " << porosity << std::endl;
std::cout << " particle_max_z: " << c_max_z << std::endl;
}
c_max_z = max_z_finder->GetValue();
DEMSim.DoDynamics(step_size);
curr_time_of_sim += step_size;
}
std::cout << "Max Z is: " << terrain_max_z << std::endl;
//DEMSim.ShowTimingStats();
//DEMSim.ClearTimingStats();
std::cout << "Exiting..." << std::endl;
return 0;
}
