Cylinder and periodic boundary conditions

I forgot to precise: this happens independently when using yadedaily (2016.06a-24-0557faf~xenial) and the latest yade version from source code (git-54c46f3)

Hi Raphael,

I am unable to test at the moment but as far as I can remember from the code periodic boundaries are not yet supported for cylinders.

Klaus

Hi Klaus,

Thank you for the answer. Can anyone confirm that ?
If it is the case, is it complicated to implement ?

Raphael

My Christmas gift below. :)

The cylinders are ok in periodic boundary conditions overall. It has been applied to sheared fiber suspensions for instance, where basically every fiber can have different nodes in different periods (if we admit that "being in different periods" makes any sense...). François also generated periodic grids IIRC.

However having bodies larger than period is the door open to remaining bugs in interaction detection, that's most likely the source of this new problem.
If you create a cylinder with length equal to period and discretized in three segments then none of the segments are too large and everything is ok. And you can still use the "largeBody" trick to insert the floor.

Bruno
__________
from yade import pack
from yade.gridpfacet import *

dp= 6e-3

periodic = 1
cylinderBottom = 0
cylindersBottom = 1

ground = 0.2*dp
H = 60*dp
L = 20*dp

O.engines = [
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_GridConnection_Aabb()],allowBiggerThanPeriod = True),
 InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Sphere_GridConnection_ScGridCoGeom()],
    [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
    [Law2_ScGeom_ViscElPhys_Basic()]
# [Ip2_FrictMat_FrictMat_FrictPhys()],
# [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 PyRunner(command = 'check()', virtPeriod = 0.1),
 NewtonIntegrator(damping=0., gravity = (0,0,-9.81))
 ]

#Material creation
O.materials.append(ViscElMat(en=0.5, et=1., young=5e6, poisson=0.5, density=2500., frictionAngle=0.4, label='Mat'))
#O.materials.append(FrictMat(young = 5e6, poisson = 0.5, density=2500,frictionAngle=0.4, label='Mat'))

#Periodic Cell or containing box
if periodic:
 O.periodic = True
 O.cell.setBox(L,L,H)
 O.bodies.append(box(center= (L/2.,L/2.,ground),extents=(100,100,0),fixed=True,color = (0.,1.,0.),wire = True,material = 'Mat'))
else:
 O.bodies.append(box(center= (L/2.,L/2.,ground),extents=(100,100,0),wire = True,fixed=True,color = (0.,1.,0.),material = 'Mat'))#Made invisible to see below
 O.bodies.append(box(center= (0,L/2.,H/2.),extents=(0,L/2.,H/2.),wire = True,fixed=True,color = (0.,1.,0.),material = 'Mat'))#Made invisible to see
 O.bodies.append(box(center= (L,L/2.,H/2.),extents=(0,L/2.,H/2.),fixed=True,color = (0.,1.,0.),material = 'Mat'))
 O.bodies.append(box(center= (L/2.,0,H/2.),extents=(L/2.,0,H/2.),wire = True,fixed=True,color = (0.,1.,0.),material = 'Mat'))#Made invisible to see
 O.bodies.append(box(center= (L/2.,L,H/2.),extents=(L/2.,0,H/2.),fixed=True,color = (0.,1.,0.),material = 'Mat'))

#Bottom fixed cylinder or fixed particles
if cylinderBottom:
 n = len(O.bodies)
 cylinder(begin = (L/2.,100*L+0.001*dp,ground+dp/2.),end =(L/2.,-100*L+0.001*dp,ground+dp/2.),radius = dp/2.,fixed = True,color = (0,0,1),intMaterial = 'Mat',extMaterial = 'Mat')#Made invisible to see inside
else:
 if cylindersBottom:
  ep=0.001*dp #a small thing
  nodes=[[L/2.,ep,ground+dp/2.],[L/2.,0.33*L+ep,ground+dp/2.],[L/2.,0.33*2*L+ep,ground+dp/2.],[L/2.,L+ep,ground+dp/2.]]
  cylinderConnection([[L/2.,ep,ground+dp/2.],[L/2.,0.33*L+ep,ground+dp/2.],[L/2.,0.33*2*L+ep,ground+dp/2.],[L/2.,L+ep,ground+dp/2.]],dp/2.,fixed = True)
 else:
  for n in range(0,int(L/dp)+1):
   O.bodies.append(sphere(center = (L/2.,n*dp,ground+dp/2.),radius = dp/2.,fixed = True,color = (0,0,1),wire = True,material = 'Mat'))#Made invisible to see inside

#Particle cloud for gravity deposition
partCloud = pack.SpherePack()
partCloud.makeCloud(minCorner=(4*L/10.,4*L/10.,ground+3.*dp),maxCorner=(6*L/10.,6*L/10.,0.9*H),rRelFuzz=0., rMean=dp/2., num = 200)
partCloud.toSimulation(material='Mat')

O.dt = 1e-5
O.saveTmp()

#Function to check if the center of a particle is contained inside the cylinder or pseudo cylinder
def check():
 for b in O.bodies:
  if b.dynamic and b.state.pos[2]<ground+dp:
   posCylXZ = Vector3(L/2.,0.,ground+dp/2.)
   pRelX = b.state.pos[0]%L - posCylXZ[0]
   pRelZ = b.state.pos[2] - posCylXZ[2]
   if sqrt(pow(pRelX,2) + pow(pRelZ,2))<dp/2.:#If the particle center is closer than the radius (i.e. if it is completely inside!)
    delta = dp/2.-sqrt(pow(pRelX,2) + pow(pRelZ,2))
    print('particle {0} is inside the cylinder from an amount of {1} diameter'.format(b.id,delta/dp))
    O.pause()

Hi Bruno,

Thanks for the answer and for the trick to split the cylinder in 3. It works very well when considering all element FrictMat and the Law2_ScGeom_FrictPhys_CundallStrack. However, when I use ViscElMat and Law2_ScGeom_ViscElPhys_Basic, I still get the same problem. Is it also the case for you ?

A pragmatic solution is to use a cylinder that is FrictMat and interact with the particles with Law2_ScGeom_FrictPhys_CundallStrack, then it works.
However, it would be nice to understand why it is not working with the Law2_ScGeom_ViscElPhys_Basic contact law, and try to fix this problem.

I will check that further but if you have any ideas in the mean time, it is more than welcome !

Raphael

Hi,

Even though I am not quite sure to understand everything, I found the origin of the problem by comparing the Law2_ScGeom_ViscElPhys_Basic to the Law2_ScGeom_FrictPhys_CundallStrack. In the following of the message, I focus only on the elastic part of Law2_ScGeom_ViscElPhys_Basic, considering a case with zero viscous damping.

When it is not periodic, there are not much difference in the evaluation of the normal and shear elastic contact force between the two contact laws, even if sign conventions are different. There is only a difference in the application of the contact force to the interacting elements.
In the Law2_ScGeom_ViscElPhys_Basic, force and torques are evaluated and directly added to the particles following:
force = phys.normalForce + shearForce + shearForceVisc;
torque1 = -c1x.cross(force);
torque2 = c2x.cross(force);

with c1x = (geom.contactPoint - de1.pos) defined between the position of the particle (de1.pos) and the contact point.
and then

addForce (id1,-force,scene);
addForce (id2, force,scene);
addTorque(id1, torque1,scene);
addTorque(id2, torque2,scene);

In the Law2_ScGeom_FrictPhys_CundallStrack it is done as:
 if (!scene->isPeriodic && !sphericalBodies) {
  State* de1 = Body::byId(id1,scene)->state.get();
  State* de2 = Body::byId(id2,scene)->state.get();
  applyForceAtContactPoint(-phys->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);}
 else {//we need to use correct branches in the periodic case, the following apply for spheres only
  Vector3r force = -phys->normalForce-shearForce;
  scene->forces.addForce(id1,force);
  scene->forces.addForce(id2,-force);
  scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
  scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
 }

Different things regarding that:
- I didn't succeed to find the location of the function applyForceAtContactPoint, but I believe that it is doing exactly the same as in the viscous law where the force and torque are applied taking into account the branch vectors c1x and c2x defined with respect to the contact point.
- We can see that when the scene is periodic (or when we consider spheres), the force is however not applied with respect to the contact point in the Law2_ScGeom_FrictPhys_CundallStrack, but with respect to the radius and the penetration depth. I think these two methods are completely equivalent for two spheres of the same size.

If I am right with these two points, there are no differences between the two contact laws regarding the elastic contribution in non-periodic frameworks and when considering monodisperse particles.

However, in the periodic case, while there is no change in the evaluation of the branch vector in Law2_ScGeom_ViscElPhys_Basic, it is different in the Law2_ScGeom_FrictPhys_CundallStrack where it is assumed that the second part of the if statement apply for all interactions. While this is a priori not appropriate for all interactions (such as an interaction between a sphere and a wall), this formulation seems to avoid some troubles. Indeed, with the Law2_ScGeom_ViscElPhys_Basic not modified formulation, the partices tend to spin a lot leading to unphysical behavior.

Could someone explain me why this trick is used for periodic BC and how it solves this spinning problem ?

Replacing the branch vector formulation based on the contact point by the one based on the penetration depth such as done in Law2_ScGeom_FrictPhys_CundallStrack, solves the whole problem, even when re-activating the viscous damping.

Last remark: The way the contact force is evaluated in Law2_ScGeom_ViscElPhys_Basic does not take into account the implemented ratcheting effect correction, and I believe it would be better to base it on the same formulation as the CundallStrack contact law to account for that. I can do it, but in order to do that, I need to know how to access the function getIncidentVel of ScGeom inside the contact law. Can someone tell me how to do that ?

Thanks !
Cheers,

Raphael

Hi,
In short you found that Law2_ScGeom_ViscElPhys_Basic does not implement periodicity and that it uses inexact expressions of the relative displacement.
I can confirm that. :)

The problem with the branch vector is because when two computational particles interact they have centers which can be arbitrarily far away - say, in different periodic windows of the periodic space. A consistent interaction geometry thus requires in principle to shift one of them to the same period where the other one is (so that (center - contactPoint) is a consistent branch vector). It is done for instance in [1]:
 applyForceAtContactPoint(-phys->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position + (scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero()));

However if the branch vector can be defined independently (like with just radius and normal, the case of CundallStrack) there is no need for a shift. The two approaches are mathematically equivalent.

You can access interaction->geom->getIncidentVel(...) where CundallStrack access interaction->geom->shearInc.
However shearInc is a cached value while getIncidentVel triggers complex calculations and needs input parameters (with possible mistakes); so I would suggest to use shearInc/scene->dt instead of getIncidentVel(...).

If improvements are in view, it may be much simpler and less energy-dissipative to introduce the viscous term directly into the CundallStrack family of functors (and leave ViscElPhys_Basic behind as it is). Since you only found "some" of the problems with this law functor.
Cheers
Bruno

[1] https://github.com/yade/trunk/blob/master/pkg/dem/CohesiveFrictionalContactLaw.cpp#L163

Hi Bruno,

Not sure I completely agree with you on all what you are saying.
- Regarding getIncidentVel, I thought at using shearInc but this gives me only the shear component and I also want the normal one so it does not work. I need to use getIncidentVel
- Regarding the periodicity, it is implemented in Law2_ScGeom_ViscElPhys_Basic. Indeed, the shift which is necessary to evaluate the geometry of the contact between two particles (scene->cell->intrShiftPos(cellDist), usually called shift2 in the code) and that you are describing, is implemented in it.
However, the problem does not seem to rely on the periodicity, at least for spheres, but on the periodicity for interaction between spheres and other objects. This can be very well identified by comparing the contact law with Law2_ScGeom_FrictPhys_CundallStrack. In the latter case, a "trick"/an approximation is used when considering periodic boundary conditions. It is the one I reported:
if (!scene->isPeriodic && !sphericalBodies) {
  State* de1 = Body::byId(id1,scene)->state.get();
  State* de2 = Body::byId(id2,scene)->state.get();
  applyForceAtContactPoint(-phys->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);}
 else {//we need to use correct branches in the periodic case, the following apply for spheres only
  Vector3r force = -phys->normalForce-shearForce;
  scene->forces.addForce(id1,force);
  scene->forces.addForce(id2,-force);
  scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
  scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
 }
Where we can see that in the case where the conditions are periodic the first condition is never fulfilled, and the part following the else is used. In that case, what is applied is that the contact point is always considered to be situated at the middle of the penetration depth. While this is ok for particles of the same size, this to my opinion does not apply when considering the interaction of a particle with a wall. My question was (and still is): why is it necessary to use this trick/approximation to avoid the bug ?
- For the fact to implement a potential correction (which I will probably do) in CundallStrack, I am not sure it is relevant, especially considering that to my knowledge a non-negligible number of people use this contact law and I believe it is not a good idea to leave this contact law with a mistake. I could simply adjust the formulation to make it similar to CundallStrack, in order to have a coherence and correct the mistake (if there is any real one).

Cheers
Raphael

> I also want the normal one

Very easy to find. Will you really waste cpu time in calculating the incident shear twice when you actually only need the normal one (you already have the shear one with my suggestion).

>contact point is always considered to be situated at the middle of the penetration depth. While this is ok for particles of the >same size, this to my opinion does not apply when considering the interaction of a particle with a wall. My question was (and >still is): why is it necessary to use this trick/approximation to avoid the bug ?

I'm not following you fully. 1/ I don't see a trick/approximation, to me it is just ordinary code. So when you ask "why" I tend to just think "why not?". 2/ why do you think it is irrelevant to include the 0.5*un term in the branch between e.g. sphere and box or spheres of different sizes? How would you like it? Did you realize that the same 0.5 is also used between facets and spheres whatever the law functor?

> I believe it is not a good idea to leave this contact law with a mistake

You are right there is nothing wrong in fixing existing code and you are very welcome to do it.
However, guess why ViscElPhys_Basic has a bug that CundallStrack does not have? Because they are two different functors instead of one, i.e. more lines of code with less eyes/fingers for each of them. Hence if someone wants to do something very useful for the years to come he could merge some of the existing functors. That was my point.
Cheers
Bruno

Hi Bruno,

Thanks for your answer, there are still some points that are unclear to me, I probably didn't understand everything:
- How do you get the normal component of the relative velocity only from the shear one ? I would need to evaluate the relative velocity for that and if I want to account for the ratcheting effect correction, it seems to me that I need to call again getIncidentVel, maybe I missed something.

>I'm not following you fully. 1/ I don't see a trick/approximation, to me it is just ordinary code. So when you ask "why" I tend to just think "why not?". 2/ why do you think it is irrelevant to include the 0.5*un term in the branch between e.g. sphere and box or spheres of different sizes? How would you like it? Did you realize that the same 0.5 is also used between facets and spheres whatever the law functor?

From what I understand the contact point is then always considered to be at half the penetration depth. Two things. 1/Why not from a modeling point of view as you say, but from a physical point of view I believe that the question of the contact point is not that simple and is fundamental. If you consider an interaction between two objects that have a different curvature, the definition of the contact point is not straightforward and you should consider the position of the contact point which gives you the right physical behavior for the normal and tangential contact force. Defining it as half the penetration depth is a choice (which I do not know if it is well reproducing the "physical"/expected behavior). 2/ If you say that the contact point is always considered as half the penetration depth, then I do not understand why that:
const State& de1 = *static_cast<State*>(bodies[id1]->state.get());
const State& de2 = *static_cast<State*>(bodies[id2]->state.get());
const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero();
const Vector3r c1x = (geom.contactPoint - de1.pos);
const Vector3r c2x = (geom.contactPoint - de2.pos - shift2);
gives different results from that:
const Vector3r c1x = geom->radius1-0.5*geom->penetrationDepth*geom->normal;
const Vector3r c2x = -geom->radius2-0.5*geom->penetrationDepth*geom->normal;

?

> However, guess why ViscElPhys_Basic has a bug that CundallStrack does not have? Because they are two different functors instead of one, i.e. more lines of code with less eyes/fingers for each of them. Hence if someone wants to do something very useful for the years to come he could merge some of the existing functors. That was my point.

I am ok with that, but this would mean keeping all the compatibility with the existing ViscElPhys_Basic law (prescription possible of (young, poisson) or (kn, ks) with (en, et) or (cn, cs)) in the new one and suppressing the ViscElPhys_Basic contact law. Is this possible ? In which file(s) should it be included ?

Raphael

> How do you get the normal component of the relative velocity?

Something like this (it would need another term from the velocity gradient, you can get inspiration from scgeom code):
(state2->vel - state1->vel).dot(normal)

>from a physical point of view I believe that the question of the contact point is not that simple and is fundamental.

Of course it is not really half-way between non-symmetric objects (different sizes, different properties, ...) but the real question is where should it be if it is not half-way?
Fundamentally, in fact, the notion of contact point/force is in itself irrelevant, that is a reason why I don't like very much the grammar " c1x = (geom.contactPoint - de1.pos)". I prefer to write equations without contact points, even if of course it can be mathematically equivalent.

The difference between the two sets of equation is that the first one uses positions, the second doesn't.
The problem in the first one is this (I didn't realize immediatly): if one body is larger than one period it is not clear how to define I->cellDist since it is not unique. If positions are not used the problem disappear.

> but this would mean keeping all the compatibility with the existing ViscElPhys_Basic law (prescription possible of (young, poisson) or (kn, ks) with (en, et) or (cn, cs))

"Merging" A and B does not mean to break B to make it fit in A by brute force. ;)
Of course the point is to keep the good stuff, else no point merging.

I guess the straight way is to make a viscous version inheriting from CundallStrack. It may need to reorganize a little bit the code of CundallStrack to let sections of the go() function be used as standalone functions (mostly a way to avoid code duplication in the inheriting classes).
Not sure yet if the law functor is enough or if some Ig/Ip functors would need something as well.