diff --git a/examples/modular/linear_spring_dashpot.py b/examples/modular/linear_spring_dashpot.py
new file mode 100644
index 0000000000000000000000000000000000000000..342af5ef82e8fba0ffefc04e04e7b4a3981db263
--- /dev/null
+++ b/examples/modular/linear_spring_dashpot.py
@@ -0,0 +1,158 @@
+import math
+import pairs
+import sys
+
+
+def update_mass_and_inertia(i):
+ rotation_matrix[i] = diagonal_matrix(1.0)
+ rotation_quat[i] = default_quaternion()
+
+ if is_sphere(i):
+ inv_inertia[i] = inversed(diagonal_matrix(0.4 * mass[i] * radius[i] * radius[i]))
+
+ else:
+ mass[i] = infinity
+ inv_inertia[i] = 0.0
+
+
+def linear_spring_dashpot(i, j):
+ delta_ij = -penetration_depth(i, j)
+ skip_when(delta_ij < 0.0)
+
+ meff = 1.0 / ((1.0 / mass[i]) + (1.0 / mass[j]))
+ stiffness_norm = meff * (pi * pi + lnDryResCoeff * lnDryResCoeff) / \
+ (collisionTime_SI * collisionTime_SI)
+ stiffness_tan = kappa * stiffness_norm
+ damping_norm = -2.0 * meff * lnDryResCoeff / collisionTime_SI
+ damping_tan = sqrt(kappa) * damping_norm
+
+ velocity_wf_i = linear_velocity[i] + cross(angular_velocity[i], contact_point(i, j) - position[i])
+ velocity_wf_j = linear_velocity[j] + cross(angular_velocity[j], contact_point(i, j) - position[j])
+
+ rel_vel = -(velocity_wf_i - velocity_wf_j)
+ rel_vel_n = dot(rel_vel, contact_normal(i, j)) * contact_normal(i, j)
+ rel_vel_t = rel_vel - rel_vel_n
+ fN = stiffness_norm * delta_ij * contact_normal(i, j) + damping_norm * rel_vel_n;
+
+ tan_spring_disp = tangential_spring_displacement[i, j]
+ impact_vel_magnitude = impact_velocity_magnitude[i, j]
+ impact_magnitude = select(impact_vel_magnitude > 0.0, impact_vel_magnitude, length(rel_vel))
+ sticking = is_sticking[i, j]
+
+ rot_tan_disp = tan_spring_disp - contact_normal(i, j) * dot(tan_spring_disp, contact_normal(i, j))
+ rot_tan_disp_len2 = squared_length(rot_tan_disp)
+ new_tan_spring_disp = dt * rel_vel_t + \
+ select(rot_tan_disp_len2 <= 0.0,
+ zero_vector(),
+ rot_tan_disp * sqrt(squared_length(tan_spring_disp) / rot_tan_disp_len2))
+
+ fTLS = stiffness_tan * new_tan_spring_disp + damping_tan * rel_vel_t
+ fTLS_len = length(fTLS)
+ t = normalized(fTLS)
+
+ f_friction_abs_static = friction_static[i, j] * length(fN)
+ f_friction_abs_dynamic = friction_dynamic[i, j] * length(fN)
+ tan_vel_threshold = 1e-8
+
+ cond1 = sticking == 1 and length(rel_vel_t) < tan_vel_threshold and fTLS_len < f_friction_abs_static
+ cond2 = sticking == 1 and fTLS_len < f_friction_abs_dynamic
+ f_friction_abs = select(cond1, f_friction_abs_static, f_friction_abs_dynamic)
+ n_sticking = select(cond1 or cond2 or fTLS_len < f_friction_abs_dynamic, 1, 0)
+ tangential_spring_displacement[i, j] = \
+ select(not cond1 and not cond2 and stiffness_tan > 0.0,
+ (f_friction_abs * t - damping_tan * rel_vel_t) / stiffness_tan,
+ new_tan_spring_disp)
+
+ impact_velocity_magnitude[i, j] = impact_magnitude
+ is_sticking[i, j] = n_sticking
+
+ fTabs = min(fTLS_len, f_friction_abs)
+ fT = fTabs * t
+ partial_force = fN + fT
+
+ apply(force, partial_force)
+ apply(torque, cross(contact_point(i, j) - position, partial_force))
+
+
+def euler(i):
+ skip_when(is_fixed(i) or is_infinite(i))
+ inv_mass = 1.0 / mass[i]
+ position[i] += 0.5 * inv_mass * force[i] * dt * dt + linear_velocity[i] * dt
+ linear_velocity[i] += inv_mass * force[i] * dt
+ wdot = rotation_matrix[i] * (inv_inertia[i] * torque[i]) * transposed(rotation_matrix[i])
+ phi = angular_velocity[i] * dt + 0.5 * wdot * dt * dt
+ rotation_quat[i] = quaternion(phi, length(phi)) * rotation_quat[i]
+ rotation_matrix[i] = quaternion_to_rotation_matrix(rotation_quat[i])
+ angular_velocity[i] += wdot * dt
+
+
+def gravity(i):
+ force[i][2] -= mass[i] * gravity_SI
+
+
+cmd = sys.argv[0]
+target = sys.argv[1] if len(sys.argv[1]) > 1 else "none"
+if target != 'cpu' and target != 'gpu':
+ print(f"Invalid target, use {cmd} <cpu/gpu>")
+
+# Config file parameters
+gravity_SI = 9.81
+ntypes = 1
+frictionCoefficient = 0.5
+restitutionCoefficient = 0.1
+poissonsRatio = 0.22
+kappa = 2.0 * (1.0 - poissonsRatio) / (2.0 - poissonsRatio) # from Thornton et al
+lnDryResCoeff = math.log(restitutionCoefficient)
+frictionStatic = 0.0
+frictionDynamic = frictionCoefficient
+
+psim = pairs.simulation(
+ "linear_spring_dashpot",
+ [pairs.sphere(), pairs.halfspace()],
+ double_prec=True,
+ use_contact_history=True,
+ particle_capacity=1000000,
+ neighbor_capacity=20)
+
+if target == 'gpu':
+ psim.target(pairs.target_gpu())
+else:
+ psim.target(pairs.target_cpu())
+ #psim.target(pairs.target_cpu(parallel=True))
+
+psim.add_position('position')
+psim.add_property('mass', pairs.real())
+psim.add_property('linear_velocity', pairs.vector())
+psim.add_property('angular_velocity', pairs.vector())
+psim.add_property('force', pairs.vector(), volatile=True)
+psim.add_property('torque', pairs.vector(), volatile=True)
+psim.add_property('radius', pairs.real())
+psim.add_property('normal', pairs.vector())
+psim.add_property('inv_inertia', pairs.matrix())
+psim.add_property('rotation_matrix', pairs.matrix())
+psim.add_property('rotation_quat', pairs.quaternion())
+
+psim.add_feature('type', ntypes)
+psim.add_feature_property('type', 'friction_static', pairs.real(), [frictionStatic for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'friction_dynamic', pairs.real(), [frictionDynamic for i in range(ntypes * ntypes)])
+
+psim.add_contact_property('is_sticking', pairs.int32(), 0)
+psim.add_contact_property('tangential_spring_displacement', pairs.vector(), [0.0, 0.0, 0.0])
+psim.add_contact_property('impact_velocity_magnitude', pairs.real(), 0.0)
+
+psim.set_domain_partitioner(pairs.regular_domain_partitioner())
+psim.pbc([False, False, False])
+psim.build_cell_lists()
+
+psim.compute(update_mass_and_inertia, symbols={'infinity': math.inf })
+
+psim.compute(gravity, symbols={'gravity_SI': gravity_SI})
+
+psim.compute(linear_spring_dashpot,
+ parameters={'dt': pairs.real(),'collisionTime_SI': pairs.real()},
+ symbols={'pi': math.pi,
+ 'kappa': kappa,
+ 'lnDryResCoeff': lnDryResCoeff})
+
+psim.compute(euler, parameters={'dt': pairs.real()})
+psim.generate()
diff --git a/examples/modular/lsd_1.cpp b/examples/modular/lsd_1.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9603f29c2abb54c66f661b4d75edfcc6f79067f4
--- /dev/null
+++ b/examples/modular/lsd_1.cpp
@@ -0,0 +1,61 @@
+#include <iostream>
+#include <memory>
+
+#include "linear_spring_dashpot.hpp"
+// cmake -DINPUT_SCRIPT=../examples/modular/linear_spring_dashpot.py -DCOMPILE_CUDA=ON -DUSE_WALBERLA=OFF -DCMAKE_BUILD_TYPE=Release -DUSER_SOURCE_FILES=../examples/modular/lsd_1.cpp ..
+
+int main(int argc, char **argv) {
+
+ auto pairs_sim = std::make_shared<PairsSimulation>();
+ pairs_sim->initialize();
+
+ auto pairs_runtime = pairs_sim->getPairsRuntime();
+
+ double dom_size = 0.1;
+ pairs_runtime->initDomain(&argc, &argv, 0, 0, 0, dom_size, dom_size, dom_size);
+
+ pairs::create_halfspace(pairs_runtime, 0,0,0, 1, 0, 0, 0, 13);
+ pairs::create_halfspace(pairs_runtime, 0,0,0, 0, 1, 0, 0, 13);
+ pairs::create_halfspace(pairs_runtime, 0,0,0, 0, 0, 1, 0, 13);
+ pairs::create_halfspace(pairs_runtime, dom_size, dom_size, dom_size, -1, 0, 0, 0, 13);
+ pairs::create_halfspace(pairs_runtime, dom_size, dom_size, dom_size, 0, -1, 0, 0, 13);
+ pairs::create_halfspace(pairs_runtime, dom_size, dom_size, dom_size, 0, 0, -1, 0, 13);
+
+ double density = 2550;
+ double diameter = 0.0029;
+ double sphere_spacing = diameter * 1.2;
+ double init_vel = 0.2;
+
+ pairs::dem_sc_grid(pairs_runtime, dom_size/2.0, dom_size/2.0, dom_size/1.5,
+ sphere_spacing,
+ diameter, diameter, diameter, init_vel, density, 1);
+
+ // pairs::create_sphere(pairs_runtime, 0.03, 0.03, 0.02, -0.2, -0.2, 0, density, diameter/2.0, 0, 0);
+ // pairs::create_sphere(pairs_runtime, 0.02, 0.02, 0.02, 0.2, 0.2, 0, density, diameter/2.0, 0, 0);
+
+ pairs_sim->update_mass_and_inertia();
+
+ pairs_sim->setCellWidth(diameter, diameter, diameter);
+ pairs_sim->setInteractionRadius(diameter);
+ pairs_sim->updateDomain();
+
+ int num_timesteps = 20000;
+ int vtk_freq = 200;
+ double dt = 1e-5;
+ double collision_time = dt*10;
+ pairs::vtk_write_subdom(pairs_runtime, "output/subdom_init", 0);
+
+ for (int t=0; t<num_timesteps; ++t){
+ if ((t%500==0) && pairs_sim->rank()==0) std::cout << "Timestep: " << t << std::endl;
+
+ pairs_sim->gravity();
+ pairs_sim->linear_spring_dashpot(collision_time, dt);
+ pairs_sim->euler(dt);
+ pairs_sim->reneighbor();
+
+ pairs::vtk_write_data(pairs_runtime, "output/sd_1_local", 0, pairs_sim->nlocal(), t, vtk_freq);
+ pairs::vtk_write_data(pairs_runtime, "output/sd_1_ghost", pairs_sim->nlocal(), pairs_sim->size(), t, vtk_freq);
+ }
+
+ pairs_sim->end();
+}
diff --git a/runtime/domain/regular_6d_stencil.cpp b/runtime/domain/regular_6d_stencil.cpp
index b469af0170a6437c4ee6c46278ebaa43ffabf4a0..d1d12f6e1dc82de3dac265b331f9b9c8cc27a90d 100644
--- a/runtime/domain/regular_6d_stencil.cpp
+++ b/runtime/domain/regular_6d_stencil.cpp
@@ -89,8 +89,8 @@ void Regular6DStencil::initialize(int *argc, char ***argv) {
this->setBoundingBox();
if (rank==0) {
std::cout << "Domain Partitioner: Regular-6D" << std::endl;
- std::cout << "Domain: [ <" << subdom_min[0] << "," << subdom_min[1] << "," << subdom_min[2] << ">, <"
- << subdom_max[0] << "," << subdom_max[1] << "," << subdom_max[2] << "> ]"<< std::endl;
+ std::cout << "Domain: [ <" << grid_min[0] << "," << grid_min[1] << "," << grid_min[2] << ">, <"
+ << grid_max[0] << "," << grid_max[1] << "," << grid_max[2] << "> ]"<< std::endl;
std::cout << "Configuration: <" << nranks[0] << "," << nranks[1] << "," << nranks[2] << ">" <<std::endl;
}
}
diff --git a/src/pairs/code_gen/interface.py b/src/pairs/code_gen/interface.py
index 9a67d659c1734c7b15bce67875023cb016546f28..bbc1699f5df5917bedc0c1934170b1f41e2775d1 100644
--- a/src/pairs/code_gen/interface.py
+++ b/src/pairs/code_gen/interface.py
@@ -34,12 +34,6 @@ class InterfaceModules:
self.refreshGhosts()
self.reverseCommunicate()
self.resetVolatiles()
-
- if self.sim._use_contact_history:
- if self.neighbor_lists:
- self.buildContactHistory(self.sim.reneighbor_frequency)
- self.resetContactHistory()
-
self.rank()
self.nlocal()
self.nghost()
@@ -126,7 +120,8 @@ class InterfaceModules:
@pairs_interface_block
def reneighbor(self):
self.sim.module_name('reneighbor')
- reneighboring_procedures = Block.from_list(self.sim, [
+
+ reneighboring_procedures = [
Exchange(self.sim._comm),
Borders(self.sim._comm),
# Note: DomainUpdateLocal must happen after exchange since local particles must be contained in AABBs
@@ -134,8 +129,16 @@ class InterfaceModules:
BuildCellListsStencil(self.sim, self.sim.cell_lists),
self.sim.update_cells_procedures,
ResetVolatileProperties(self.sim)
- ])
- self.sim.add_statement(reneighboring_procedures)
+ ]
+
+ if self.sim._use_contact_history:
+ reneighboring_procedures += [
+ BuildContactHistory(self.sim, self.sim._contact_history, self.sim.cell_lists),
+ ResetContactHistoryUsageStatus(self.sim, self.sim._contact_history),
+ ClearUnusedContactHistory(self.sim, self.sim._contact_history)
+ ]
+
+ self.sim.add_statement(Block.from_list(self.sim, reneighboring_procedures))
@pairs_interface_block
def refreshGhosts(self):
@@ -152,25 +155,6 @@ class InterfaceModules:
self.sim.module_name('resetVolatiles')
self.sim.add_statement(ResetVolatileProperties(self.sim))
- @pairs_interface_block
- def buildContactHistory(self, reneighbor_frequency=1):
- self.sim.module_name('buildContactHistory')
- timestep = Parameter(self.sim, f'timestep', Types.Int32)
- cond = ScalarOp.inline(ScalarOp.or_op(
- ScalarOp.cmp((timestep + 1) % reneighbor_frequency, 0),
- ScalarOp.cmp(timestep, 0)
- ))
-
- self.sim.add_statement(
- Filter(self.sim, cond,
- BuildContactHistory(self.sim, self.sim._contact_history, self.sim.cell_lists)))
-
- @pairs_interface_block
- def resetContactHistory(self):
- self.sim.module_name('resetContactHistory')
- self.sim.add_statement(ResetContactHistoryUsageStatus(self.sim, self.sim._contact_history))
- self.sim.add_statement(ClearUnusedContactHistory(self.sim, self.sim._contact_history))
-
@pairs_interface_block
def rank(self):
self.sim.module_name('rank')
diff --git a/src/pairs/sim/interaction.py b/src/pairs/sim/interaction.py
index 2671a907e526e1b9aefe0677349c5d2233f962f7..49aff87d50c1f6b7d508e54b3df3347d1004869d 100644
--- a/src/pairs/sim/interaction.py
+++ b/src/pairs/sim/interaction.py
@@ -47,8 +47,8 @@ class NeighborFor:
self.particle = particle
self.cell_lists = cell_lists
self.neighbor_lists = neighbor_lists
- # self.shapes = range(sim.max_shapes()) if shapes is None else shapes
- self.shapes = [shapes]
+ self.shapes = range(sim.max_shapes()) if shapes is None else [shapes]
+ # self.shapes = [shapes]
def __str__(self):
return f"NeighborFor<{self.particle}>"