diff --git a/examples/dem_sd.py b/examples/dem_sd.py
new file mode 100644
index 0000000000000000000000000000000000000000..737da4b562249d089628bdf12926d833396f7415
--- /dev/null
+++ b/examples/dem_sd.py
@@ -0,0 +1,177 @@
+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)
+ 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))
+ rel_vel_t = rel_vel - rel_vel_n * contact_normal(i, j)
+ fNabs = stiffness_norm * delta_ij + damping_norm * rel_vel_n
+ fN = fNabs * contact_normal(i, j)
+
+ fTabs = min(damping_tan * length(rel_vel_t), friction_dynamic[i, j] * fNabs) # static or dynamic?
+ fT = fTabs * normalized(rel_vel_t)
+
+ partial_force = fN + fT
+
+ apply(force, partial_force)
+ apply(torque, cross(contact_point(i, j) - position, partial_force))
+
+
+def euler(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):
+ volume = (4.0 / 3.0) * pi * radius[i] * radius[i] * radius[i]
+ force[i][2] = force[i][2] - (densityParticle_SI - densityFluid_SI) * volume * gravity_SI
+ # if type[i]==0:
+ # printf("gravity", i, volume, force[i][2], volume/2, 2>1)
+
+
+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
+domainSize_SI=[0.1, 0.1, 0.1]
+diameter_SI = 0.009
+gravity_SI = 9.81
+densityFluid_SI = 1000
+densityParticle_SI = 2550
+generationSpacing_SI = 0.01
+initialVelocity_SI = 1
+dt_SI = 5e-5
+frictionCoefficient = 0.5
+restitutionCoefficient = 0.1
+collisionTime_SI = 5e-4
+poissonsRatio = 0.22
+timeSteps = 10000
+visSpacing = 200
+denseBottomLayer = False
+bottomLayerOffsetFactor = 1.0
+kappa = 2.0 * (1.0 - poissonsRatio) / (2.0 - poissonsRatio) # from Thornton et al
+
+minDiameter_SI = diameter_SI #* 0.9
+maxDiameter_SI = diameter_SI #* 1.1
+linkedCellWidth = 1.01 * maxDiameter_SI
+
+ntypes = 1
+lnDryResCoeff = math.log(restitutionCoefficient)
+frictionStatic = 0.0
+frictionDynamic = frictionCoefficient
+
+psim = pairs.simulation(
+ "dem_sd",
+ [pairs.sphere(), pairs.halfspace()],
+ timesteps=timeSteps,
+ double_prec=True,
+ use_contact_history=False,
+ particle_capacity=1000000,
+ neighbor_capacity=20,
+ debug=True, generate_whole_program=False)
+
+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(), 1.0)
+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(), 1.0)
+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.set_domain([0.0, 0.0, 0.0, domainSize_SI[0], domainSize_SI[1], domainSize_SI[2]])
+# psim.set_domain_partitioner(pairs.block_forest(), initDomainFromWalberla=True)
+psim.set_domain_partitioner(pairs.block_forest())
+# psim.set_domain_partitioner(pairs.regular_domain_partitioner_xy())
+psim.pbc([True, True, False])
+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",
+# ['type', 'mass', 'radius', 'position', 'linear_velocity', 'flags'],
+# pairs.sphere())
+
+psim.read_particle_data(
+ "data/planes.input",
+ ['uid', 'type', 'mass', 'position', 'normal', 'flags'],
+ pairs.halfspace())
+
+psim.setup(update_mass_and_inertia, {'densityParticle_SI': densityParticle_SI,
+ 'pi': math.pi,
+ 'infinity': math.inf })
+
+#psim.compute_half()
+psim.build_cell_lists(linkedCellWidth)
+psim.vtk_output(f"output/dem_{target}", frequency=visSpacing)
+
+psim.compute(gravity,
+ symbols={'densityParticle_SI': densityParticle_SI,
+ 'densityFluid_SI': densityFluid_SI,
+ 'gravity_SI': gravity_SI,
+ 'pi': math.pi })
+
+psim.compute(linear_spring_dashpot,
+ linkedCellWidth,
+ symbols={'dt': dt_SI,
+ 'pi': math.pi,
+ 'kappa': kappa,
+ 'lnDryResCoeff': lnDryResCoeff,
+ 'collisionTime_SI': collisionTime_SI})
+
+psim.compute(euler, symbols={'dt': dt_SI})
+psim.generate()
diff --git a/src/pairs/ir/print.py b/src/pairs/ir/print.py
new file mode 100644
index 0000000000000000000000000000000000000000..58c2a6e1d062342e6e2d6511fccc8ee75a7ce55a
--- /dev/null
+++ b/src/pairs/ir/print.py
@@ -0,0 +1,24 @@
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.lit import Lit
+
+class Print(ASTNode):
+ def __init__(self, sim, *args):
+ super().__init__(sim)
+ self.args = [Lit.cvt(sim, a) for a in args]
+ self.sim.add_statement(self)
+
+ def children(self):
+ return self.args
+
+ def __str__(self):
+ return "Print<" + ", ".join(str(arg) for arg in self.args) + ">"
+
+class PrintCode(ASTNode):
+ def __init__(self, sim, str):
+ super().__init__(sim)
+ self.arg = Lit.cvt(sim, str)
+ self.sim.add_statement(self)
+
+ def children(self):
+ return self.arg
+
\ No newline at end of file