diff --git a/examples/dem.py b/examples/dem.py
index 2c9dcf3c5bac4294e832f5a85acaec2499ba8b01..6f426995f0b6f4607ad6eace5d99e9ecfa2f8944 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -8,7 +8,6 @@ def update_mass_and_inertia(i):
rotation_quat[i] = default_quaternion()
if is_sphere(i):
- #mass[i] = ((4.0 / 3.0) * pi) * radius[i] * radius[i] * radius[i] * densityParticle_SI
inv_inertia[i] = inversed(diagonal_matrix(0.4 * mass[i] * radius[i] * radius[i]))
else:
@@ -164,13 +163,17 @@ psim.add_contact_property('impact_velocity_magnitude', pairs.real(), 0.0)
psim.set_domain([0.0, 0.0, 0.0, domainSize_SI[0], domainSize_SI[1], domainSize_SI[2]])
psim.set_domain_partitioner(pairs.regular_domain_partitioner_xy())
psim.pbc([True, True, False])
-psim.read_particle_data(
- "data/spheres.input",
- #"data/spheres_4x4x2.input",
- #"data/spheres_6x6x2.input",
- #"data/spheres_8x8x2.input",
- ['uid', 'type', 'mass', 'radius', 'position', 'linear_velocity', 'flags'],
- pairs.sphere())
+psim.dem_sc_grid(
+ domainSize_SI[0], domainSize_SI[1], domainSize_SI[2], generationSpacing_SI,
+ diameter_SI, minDiameter_SI, maxDiameter_SI, initialVelocity_SI, densityParticle_SI, ntypes)
+
+#psim.read_particle_data(
+# "data/spheres.input",
+# "data/spheres_4x4x2.input",
+# "data/spheres_6x6x2.input",
+# #"data/spheres_8x8x2.input",
+# ['uid', 'type', 'mass', 'radius', 'position', 'linear_velocity', 'flags'],
+# pairs.sphere())
#psim.read_particle_data(
# "data/spheres_bottom.input",
@@ -187,8 +190,8 @@ psim.setup(update_mass_and_inertia, {'densityParticle_SI': densityParticle_SI,
'infinity': math.inf })
#psim.compute_half()
-psim.build_cell_lists(linkedCellWidth, store_neighbors_per_cell=True)
-#psim.vtk_output(f"output/dem_{target}", frequency=visSpacing)
+psim.build_cell_lists(linkedCellWidth)
+psim.vtk_output(f"output/dem_{target}", frequency=visSpacing)
psim.compute(gravity,
symbols={'densityParticle_SI': densityParticle_SI,
diff --git a/runtime/dem_sc_grid.hpp b/runtime/dem_sc_grid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6f47fc5c50cb881fd72154212f370b13649b6e24
--- /dev/null
+++ b/runtime/dem_sc_grid.hpp
@@ -0,0 +1,155 @@
+#include <iostream>
+#include <math.h>
+#include <random>
+//---
+#include "pairs.hpp"
+
+#pragma once
+
+namespace pairs {
+
+namespace internal {
+
+static std::mt19937 generator; // static std::mt19937_64 generator;
+
+std::mt19937 & get_generator() {
+ // std::mt19937_64
+ return generator;
+}
+
+}
+
+template< typename REAL_TYPE = real_t>
+REAL_TYPE realRandom(
+ const REAL_TYPE min = REAL_TYPE(0),
+ const REAL_TYPE max = REAL_TYPE(1),
+ std::mt19937& generator = internal::get_generator()) {
+
+ static_assert(
+ std::numeric_limits<REAL_TYPE>::is_specialized &&
+ !std::numeric_limits<REAL_TYPE>::is_integer,
+ "Floating point type required/expected!" );
+
+ std::uniform_real_distribution<REAL_TYPE> distribution(min, max);
+
+ REAL_TYPE value;
+#ifdef _OPENMP
+ #pragma omp critical (Random_random)
+#endif
+ { value = distribution( generator ); }
+
+ return value;
+}
+
+
+
+template<typename REAL_TYPE> class RealRandom {
+public:
+ RealRandom(const std::mt19937::result_type& seed = std::mt19937::result_type()) { generator_.seed(seed); }
+ REAL_TYPE operator()(const REAL_TYPE min = REAL_TYPE(0), const REAL_TYPE max = REAL_TYPE(1) ) {
+ return realRandom(min, max, generator_);
+ }
+private:
+ std::mt19937 generator_;
+};
+
+bool point_within_aabb(double point[], double aabb[]) {
+ return point[0] >= aabb[0] && point[0] < aabb[3] &&
+ point[1] >= aabb[1] && point[1] < aabb[4] &&
+ point[2] >= aabb[2] && point[2] < aabb[5];
+}
+
+int dem_sc_grid(PairsSimulation *ps, double xmax, double ymax, double zmax, double spacing, double diameter, double min_diameter, double max_diameter, double initial_velocity, double particle_density, int ntypes) {
+ auto uid = ps->getAsIntegerProperty(ps->getPropertyByName("uid"));
+ auto shape = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
+ auto types = ps->getAsIntegerProperty(ps->getPropertyByName("type"));
+ auto flags = ps->getAsIntegerProperty(ps->getPropertyByName("flags"));
+ auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
+ auto radius = ps->getAsFloatProperty(ps->getPropertyByName("radius"));
+ auto positions = ps->getAsVectorProperty(ps->getPropertyByName("position"));
+ auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
+ int last_uid = 1;
+ int nparticles = 0;
+
+ const double xmin = 0.0;
+ const double ymin = 0.0;
+ const double zmin = diameter;
+
+ double gen_domain[] = {xmin, ymin, zmin, xmax, ymax, zmax};
+ double ref_point[] = {spacing * 0.5, spacing * 0.5, spacing * 0.5};
+ double sc_xmin = xmin - ref_point[0];
+ double sc_ymin = ymin - ref_point[1];
+ double sc_zmin = zmin - ref_point[2];
+
+ int iret = (int)(ceil(sc_xmin / spacing));
+ int jret = (int)(ceil(sc_ymin / spacing));
+ int kret = (int)(ceil(sc_zmin / spacing));
+
+ int i = iret;
+ int j = jret;
+ int k = kret;
+
+ double point[3];
+ point[0] = ref_point[0] + i * spacing;
+ point[1] = ref_point[1] + j * spacing;
+ point[2] = ref_point[2] + k * spacing;
+
+ while(point_within_aabb(point, gen_domain)) {
+ int particle_uid = last_uid;
+ auto diameter = realRandom<real_t>(min_diameter, max_diameter);
+
+ if(ps->getDomainPartitioner()->isWithinSubdomain(point[0], point[1], point[2])) {
+ real_t rad = diameter * 0.5;
+ uid(nparticles) = particle_uid;
+ radius(nparticles) = rad;
+ masses(nparticles) = ((4.0 / 3.0) * M_PI) * rad * rad * rad * particle_density;
+ positions(nparticles, 0) = point[0];
+ positions(nparticles, 1) = point[1];
+ positions(nparticles, 2) = point[2];
+ velocities(nparticles, 0) = 0.1 * realRandom<real_t>(-initial_velocity, initial_velocity);
+ velocities(nparticles, 1) = 0.1 * realRandom<real_t>(-initial_velocity, initial_velocity);
+ velocities(nparticles, 2) = -initial_velocity;
+ types(nparticles) = rand() % ntypes;
+ flags(nparticles) = 0;
+ shape(nparticles) = 0; // sphere
+
+ std::cout << uid(nparticles) << "," << types(nparticles) << "," << masses(nparticles) << "," << radius(nparticles) << ","
+ << positions(nparticles, 0) << "," << positions(nparticles, 1) << "," << positions(nparticles, 2) << ","
+ << velocities(nparticles, 0) << "," << velocities(nparticles, 1) << "," << velocities(nparticles, 2) << ","
+ << flags(nparticles) << std::endl;
+
+ nparticles++;
+ }
+
+ ++i;
+ point[0] = ref_point[0] + i * spacing;
+ point[1] = ref_point[1] + j * spacing;
+ point[2] = ref_point[2] + k * spacing;
+
+ if(!point_within_aabb(point, gen_domain)) {
+ i = iret;
+ j++;
+ point[0] = ref_point[0] + i * spacing;
+ point[1] = ref_point[1] + j * spacing;
+ point[2] = ref_point[2] + k * spacing;
+
+ if(!point_within_aabb(point, gen_domain)) {
+ j = jret;
+ k++;
+ point[0] = ref_point[0] + i * spacing;
+ point[1] = ref_point[1] + j * spacing;
+ point[2] = ref_point[2] + k * spacing;
+
+ if(!point_within_aabb(point, gen_domain)) {
+ break;
+ }
+ }
+ }
+
+ last_uid++;
+ }
+
+ return nparticles;
+}
+
+}
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index 0d6d8b05c0bf75d410e835e22144e89514bd83a9..03c52993bf4dbf408a0ad0cce3cf99d6f425179b 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -74,6 +74,7 @@ class CGen:
self.print("//---")
self.print("#include \"runtime/likwid-marker.h\"")
self.print("#include \"runtime/copper_fcc_lattice.hpp\"")
+ self.print("#include \"runtime/dem_sc_grid.hpp\"")
self.print("#include \"runtime/pairs.hpp\"")
self.print("#include \"runtime/read_from_file.hpp\"")
self.print("#include \"runtime/timing.hpp\"")
diff --git a/src/pairs/sim/dem_sc_grid.py b/src/pairs/sim/dem_sc_grid.py
new file mode 100644
index 0000000000000000000000000000000000000000..c7ef027bf8308d748f139a06afc0953ac649b650
--- /dev/null
+++ b/src/pairs/sim/dem_sc_grid.py
@@ -0,0 +1,29 @@
+from pairs.ir.assign import Assign
+from pairs.ir.block import pairs_inline
+from pairs.ir.functions import Call_Int
+from pairs.sim.lowerable import Lowerable
+
+
+class DEMSCGrid(Lowerable):
+ def __init__(
+ self, sim, xmax, ymax, zmax, spacing, diameter, min_diameter, max_diameter, initial_velocity, particle_density, ntypes):
+
+ super().__init__(sim)
+ self._xmax = xmax
+ self._ymax = ymax
+ self._zmax = zmax
+ self._spacing = spacing
+ self._diameter = diameter
+ self._min_diameter = min_diameter
+ self._max_diameter = max_diameter
+ self._initial_velocity = initial_velocity
+ self._particle_density = particle_density
+ self._ntypes = ntypes
+ #sim.set_domain([0.0, 0.0, 0.0, self._xprd, self._yprd, self._zprd])
+
+ @pairs_inline
+ def lower(self):
+ Assign(self.sim, self.sim.nlocal,
+ Call_Int(self.sim, "pairs::dem_sc_grid",
+ [self._xmax, self._ymax, self._zmax, self._spacing, self._diameter, self._min_diameter, self._max_diameter,
+ self._initial_velocity, self._particle_density, self._ntypes]))
diff --git a/src/pairs/sim/interaction.py b/src/pairs/sim/interaction.py
index 4d3baa6d325faf614c781fd3e4f1d959a8f857ad..5502d8514e08d7f88e5d4054d2b7457a6928c7c0 100644
--- a/src/pairs/sim/interaction.py
+++ b/src/pairs/sim/interaction.py
@@ -75,8 +75,11 @@ class NeighborFor:
if self.sim._compute_half:
shape_id = self.sim.get_shape_id(shape)
- shape_cond = particle_shape[particle_id] > shape_id
- condition = ScalarOp.or_op(shape_cond, self.particle < particle_id)
+ condition = ScalarOp.or_op(
+ particle_shape[particle_id] > shape_id,
+ ScalarOp.and_op(
+ ScalarOp.cmp(particle_shape[particle_id], shape_id),
+ self.particle < particle_id))
else:
condition = ScalarOp.neq(particle_id, self.particle)
@@ -103,8 +106,11 @@ class NeighborFor:
if self.sim._compute_half:
shape_id = self.sim.get_shape_id(shape)
- shape_cond = particle_shape[particle_id] > shape_id
- condition = ScalarOp.or_op(shape_cond, self.particle < particle_id)
+ condition = ScalarOp.or_op(
+ particle_shape[particle_id] > shape_id,
+ ScalarOp.and_op(
+ ScalarOp.cmp(particle_shape[particle_id], shape_id),
+ self.particle < particle_id))
else:
condition = ScalarOp.neq(particle_id, self.particle)
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 8b10fe053220f2d7240e8b6898a2750cd16606d6..291f085666fdb6cdac4e4f3d3b3a28969ab76c78 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -16,6 +16,7 @@ from pairs.sim.cell_lists import CellLists, BuildCellLists, BuildCellListsStenci
from pairs.sim.comm import Comm
from pairs.sim.contact_history import ContactHistory, BuildContactHistory, ClearUnusedContactHistory, ResetContactHistoryUsageStatus
from pairs.sim.copper_fcc_lattice import CopperFCCLattice
+from pairs.sim.dem_sc_grid import DEMSCGrid
from pairs.sim.domain import InitializeDomain
from pairs.sim.domain_partitioners import DomainPartitioners
from pairs.sim.domain_partitioning import DimensionRanges
@@ -241,6 +242,11 @@ class Simulation:
def copper_fcc_lattice(self, nx, ny, nz, rho, temperature, ntypes):
self.setups.add_statement(CopperFCCLattice(self, nx, ny, nz, rho, temperature, ntypes))
+ def dem_sc_grid(self, xmax, ymax, zmax, spacing, diameter, min_diameter, max_diameter, initial_velocity, particle_density, ntypes):
+ self.setups.add_statement(
+ DEMSCGrid(self, xmax, ymax, zmax, spacing, diameter, min_diameter, max_diameter,
+ initial_velocity, particle_density, ntypes))
+
def build_cell_lists(self, spacing, store_neighbors_per_cell=False):
self._store_neighbors_per_cell = store_neighbors_per_cell
self.cell_lists = CellLists(self, self._dom_part, spacing, spacing)