From af61e478292f1744663efe77d4d9050a9ec33bba Mon Sep 17 00:00:00 2001 From: shouyu-korokke Date: Fri, 1 Dec 2023 23:42:31 +0100 Subject: [PATCH] Create pyDEME_Shake.py --- src/demo/pyDEME_Shake.py | 167 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100644 src/demo/pyDEME_Shake.py diff --git a/src/demo/pyDEME_Shake.py b/src/demo/pyDEME_Shake.py new file mode 100644 index 0000000..d2316a2 --- /dev/null +++ b/src/demo/pyDEME_Shake.py @@ -0,0 +1,167 @@ +# Copyright (c) 2021, SBEL GPU Development Team +# Copyright (c) 2021, University of Wisconsin - Madison +# +# SPDX-License-Identifier: BSD-3-Clause + +# ============================================================================= +# Put particles in a jar, the shake the jar in the hope to change the bulk density. +# ============================================================================= + +import DEME + +import numpy as np +import os + + +if __name__ == "__main__": + out_dir = "DemoOutput_Shake/" + out_dir = os.path.join(os.getcwd(), out_dir) + os.makedirs(out_dir, exist_ok=True) + + DEMSim = DEME.DEMSolver() + DEMSim.SetVerbosity("INFO") + DEMSim.SetOutputFormat("CSV") + DEMSim.SetOutputContent(["ABSV"]) + DEMSim.SetMeshOutputFormat("VTK") + DEMSim.SetContactOutputContent(["OWNER", "FORCE", "POINT"]) + + # If you don't need individual force information, then this option makes the solver run a bit faster. + DEMSim.SetNoForceRecord() + + # E, nu, CoR, mu, Crr... + mat_type_cone = DEMSim.LoadMaterial( + {"E": 5e7, "nu": 0.3, "CoR": 0.5}) + mat_type_terrain = DEMSim.LoadMaterial( + {"E": 5e7, "nu": 0.3, "CoR": 0.5}) + DEMSim.UseFrictionlessHertzianModel() + + + shake_amp = 0.1 + shake_speed = 2 # Num of periods per second + step_size = 1e-5 + world_size = 2. + soil_bin_diameter = 0.584 + cone_surf_area = 323e-6 + cone_diameter = np.sqrt(cone_surf_area / np.pi) * 2 + 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_terrain) + # Now add a cylinderical boundary along with a bottom plane + bottom = -0.5 + walls = DEMSim.AddExternalObject() + walls.AddCylinder([0, 0, 0], [0, 0, 1], soil_bin_diameter / 2., mat_type_terrain, 0) + walls.AddPlane([0, 0, bottom], [0, 0, 1], mat_type_terrain) + walls.SetFamily(1) + + # Define the terrain particle templates. Two types of clumps. + terrain_density = 2.6e3 + clump_vol = 5.5886717 + mass = terrain_density * clump_vol + MOI1 = np.array([2.928, 2.6029, 3.9908]) * terrain_density + clump_vol2 = 2.1670011 + mass2 = terrain_density * clump_vol2 + MOI2 = np.array([0.57402126, 0.60616378, 0.92890173]) * terrain_density + # Then load it to system + template_1 = DEMSim.LoadClumpType(mass, MOI1.tolist( + ), DEME.GetDEMEDataFile("clumps/triangular_flat.csv"), mat_type_terrain) + template_2 = DEMSim.LoadClumpType(mass, MOI1.tolist( + ), DEME.GetDEMEDataFile("clumps/triangular_flat_6comp.csv"), mat_type_terrain) + # Decide the scalings of the templates we just created (so that they are... like particles, not rocks) + scale1 = 0.015 + scale2 = 0.004 + template_1.Scale(scale1) + template_2.Scale(scale2) + + # Sampler to sample + sampler1 = DEME.GridSampler(scale1 * 3.) + fill_height = 0.5 + fill_center = [0, 0, bottom + fill_height / 2] + fill_radius = soil_bin_diameter / 2. - scale1 * 2. + input_xyz = sampler1.SampleCylinderZ( + fill_center, fill_radius, fill_height / 2 - scale1 * 2.) + DEMSim.AddClumps(template_1, input_xyz) + # Another batch... + sampler2 = DEME.GridSampler(scale2 * 3.) + fill_center += [0, 0, fill_height] + fill_radius = soil_bin_diameter / 2. - scale2 * 2. + input_xyz = sampler2.SampleCylinderZ( + fill_center, fill_radius, fill_height / 2 - scale2 * 2.) + DEMSim.AddClumps(template_2, input_xyz) + + # Now add a `cap' to the container when we shake it + compressor = DEMSim.AddExternalObject() + compressor.AddPlane([0, 0, 0], [0, 0, -1], mat_type_terrain) + compressor.SetFamily(1) + compressor_tracker = DEMSim.Track(compressor) + + # Family 1 shakes, family 2 is fixed + shake_pattern_x = DEME.to_string_with_precision(shake_amp) + " * sin(" + DEME.to_string_with_precision(shake_speed) + " * 2 * deme::PI * t)" + DEMSim.SetFamilyPrescribedMotion(1, shake_pattern_x, "0", shake_pattern_x) + DEMSim.SetFamilyFixed(2) + + # Some inspectors + max_z_finder = DEMSim.CreateInspector("clump_max_z") + min_z_finder = DEMSim.CreateInspector("clump_min_z") + total_mass_finder = DEMSim.CreateInspector("clump_mass") + max_v_finder = DEMSim.CreateInspector("clump_max_absv") + + DEMSim.SetInitTimeStep(step_size) + DEMSim.SetGravitationalAcceleration([0, 0, -9.81]) + + # 256 or 512 are common choices. Note that in cases where the force model is modified, too many registers may be + # used in the kernel, so we have to reduce this number to use 256. In other cases (and most cases), 512 is fine and + # may make the code run a bit faster. Usually, the user do not have to call SetForceCalcThreadsPerBlock if they + # don't know the implication. + DEMSim.SetForceCalcThreadsPerBlock(512) + DEMSim.Initialize() + + # Settle phase + currframe = 0 + curr_step = 0 + fps = 10 + out_steps = (int)(1.0 / (fps * step_size)) + compressor_tracker.SetPos([0, 0, max_z_finder.GetValue()]) + for i in range (0, 5): + filename = os.path.join(out_dir, f"DEMdemo_output_{currframe:04d}.csv") + DEMSim.WriteSphereFile(filename) + DEMSim.DoDynamicsThenSync(0.1) + currframe += 1 + + stop_time = 2.0 # Time before stopping shaking and measure bulk density + stop_steps = (int)(stop_time * (1.0 / step_size)) + sim_end = 6.0 + print(f"Output at {fps} FPS") + + t = 0. + while(t < sim_end): + if (curr_step % out_steps == 0): + print(f"Frame: {currframe}", flush=True) + filename = os.path.join( + out_dir, f"DEMdemo_output_{currframe:04d}.csv") + DEMSim.WriteSphereFile(filename) + print(f"Max system velocity: {max_z_finder.GetValue()}") + currframe += 1 + DEMSim.ShowThreadCollaborationStats() + if (curr_step % stop_steps == 0): + # Measure + max_z = max_z_finder.GetValue() + min_z = min_z_finder.GetValue() + matter_mass = total_mass_finder.GetValue() + total_volume = np.pi * (soil_bin_diameter * soil_bin_diameter /4) + bulk_density = matter_mass / total_volume + print(f"Max z: {max_z}") + print(f"Min z: {min_z}") + print(f"Bulk density: {bulk_density}") + # Put the cap to its new position and start shaking + compressor_tracker.SetPos([0, 0, max_z]) + + DEMSim.DoDynamics(step_size) + t += step_size + curr_step += 1 + + # Output the final configuration of the clumps as a file. This is just a demonstration. This particular + # configuration is not that useful as no other demos actually use it, unlike the GRC-1 soil. + cp_filename = os.path.join(out_dir, f"GRC_rho{bulk_density}.csv") + DEMSim.WriteSphereFile(cp_filename) + + print("Shake demo exiting...")