yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05936
[Branch ~yade-dev/yade/trunk] Rev 2501: 1. Fix kinetic energy in PBC
------------------------------------------------------------
revno: 2501
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Tue 2010-10-19 16:57:27 +0200
message:
1. Fix kinetic energy in PBC
2. Add tests to check PBC: Ek, incident velocity with ScGeom (with and without ratcheting), homothetic resize
3. Make it possible to set O.dt=0, for reasons of testing.
renamed:
py/tests/cell.py => py/tests/pbc.py
modified:
pkg/dem/Shop.cpp
py/SConscript
py/tests/__init__.py
py/tests/omega.py
py/tests/wrapper.py
py/wrapper/yadeWrapper.cpp
py/tests/pbc.py
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp 2010-10-18 19:27:07 +0000
+++ pkg/dem/Shop.cpp 2010-10-19 14:57:27 +0000
@@ -189,38 +189,10 @@
if(!b || !b->isDynamic()) continue;
const State* state(b->state.get());
// ½(mv²+ωIω)
- Vector3r vel;
- if(!isPeriodic) vel=b->state->vel;
- else{
- /* Kinetic energy is defined with regards to absolute space (or off by a constant for system in steady motion).
- However, in the periodic case with non-zero velocity gradient, there is no such inertial system,
- as velocity depends on absolute coordinate in space (NewtonIntegrator::homotheticCellResize).
- There are several inconsistent options (and no consistent one):
-
- 1. Subtract the contribution of homotheticCellResize; such kinetic energy however materializes
- when there is an interaction between particles, and it would come apparently from nowhere.
-
- 2. Compute the kinetic energy for particles translated inside the periodic cell. That way,
- we avoid the dependency on absolute space position, while still keeping the contribution
- of homothetic resize. This comes at the cost of discontinuity which is apparent for 1 particle.
-
- Suppose that particle A is inside the base cell (i.e. between the origin and O.cell.size) and that
- it has some velocity away from the origin. Applying velocity gradient, velocity of the particle will
- be augmented by the amount corrseponding to its distance from the origin. At some moment, A will
- leave the base cell. Kinetic energy computed from coordinates wrapped inside the base cell will suddenly
- drop, since the wrapped coordinate will jump from O.cell.size[i] to 0.
-
- With increasing number of particles, this effect will become less and less apparent,
- eventually approaching zero. Moreoved, for dense packings, two directions of motion accross the cell
- boundary will compensate each other.
-
- 3. Propose some special form of potential energy related to the periodic cell. Not sure if that is doable,
- since the amount of energy depends also on what interaction (with which particles) will be established
- in the future.
-
- */
- Vector3i period; Vector3r basePos=scene->cell->wrapShearedPt(state->pos,period);
- vel=state->vel-scene->cell->velGrad*(state->pos-basePos);
+ Vector3r vel=b->state->vel;
+ if(isPeriodic){
+ /* Only take in account the fluctuation velocity, not the mean velocity of homothetic resize. */
+ vel-=scene->cell->velGrad*state->pos;
// TODO: move NewtonIntegrator.homotheticCellResize to Cell.homotheticDeformation so that
// we have access to how is the velocity adjusted
// in addition, create function in Cell that will compute velocity compensations etc for us
=== modified file 'py/SConscript'
--- py/SConscript 2010-10-18 19:27:07 +0000
+++ py/SConscript 2010-10-19 14:57:27 +0000
@@ -46,7 +46,7 @@
env.File('__init__.py','tests'),
env.File('wrapper.py','tests'),
env.File('omega.py','tests'),
- env.File('cell.py','tests'),
+ env.File('pbc.py','tests'),
])
# 3rd party modules:
=== modified file 'py/tests/__init__.py'
--- py/tests/__init__.py 2010-10-18 19:27:07 +0000
+++ py/tests/__init__.py 2010-10-19 14:57:27 +0000
@@ -4,7 +4,7 @@
import unittest,inspect
# add any new test suites to the list here, so that they are picked up by testAll
-allTests=['wrapper','omega','cell']
+allTests=['wrapper','omega','pbc']
# all yade modules (ugly...)
import yade.eudoxos,yade.export,yade.linterpolation,yade.log,yade.pack,yade.plot,yade.post2d,yade.timing,yade.utils,yade.ymport
=== modified file 'py/tests/omega.py'
--- py/tests/omega.py 2010-10-18 19:27:07 +0000
+++ py/tests/omega.py 2010-10-19 14:57:27 +0000
@@ -24,28 +24,28 @@
def setUp(self):
O.reset(); O.periodic=True
def testAttributesAreCrossUpdated(self):
- "updates Hsize automatically when refSize is updated"
+ "Cell: updates Hsize automatically when refSize is updated"
O.cell.refSize=(2.55,11,45)
self.assert_(O.cell.Hsize==Matrix3(2.55,0,0, 0,11,0, 0,0,45));
class TestMaterialStateAssociativity(unittest.TestCase):
def setUp(self): O.reset()
def testThrowsAtBadCombination(self):
- "throws when body has material and state that don't work together."
+ "Material+State: throws when body has material and state that don't work together."
b=Body()
b.mat=CpmMat()
b.state=State() #should be CpmState()
O.bodies.append(b)
self.assertRaises(RuntimeError,lambda: O.step()) # throws runtime_error
def testThrowsAtNullState(self):
- "throws when body has material but NULL state."
+ "Material+State: throws when body has material but NULL state."
b=Body()
b.mat=Material()
b.state=None # → shared_ptr<State>() by boost::python
O.bodies.append(b)
self.assertRaises(RuntimeError,lambda: O.step())
def testMaterialReturnsState(self):
- "CpmMat returns CpmState when asked for newAssocState"
+ "Material+State: CpmMat returns CpmState when asked for newAssocState"
self.assert_(CpmMat().newAssocState().__class__==CpmState)
class TestBodies(unittest.TestCase):
@@ -55,22 +55,22 @@
O.bodies.append([utils.sphere([random.random(),random.random(),random.random()],random.random()) for i in range(0,self.count)])
random.seed()
def testIterate(self):
- "Iteration over O.bodies"
+ "Bodies: Iteration"
counted=0
for b in O.bodies: counted+=1
self.assert_(counted==self.count)
def testLen(self):
- "len(O.bodies)"
+ "Bodies: len(O.bodies)"
self.assert_(len(O.bodies)==self.count)
def testErase(self):
- "Erased bodies are None in python"
+ "Bodies: erased bodies are None in python"
O.bodies.erase(0)
self.assert_(O.bodies[0]==None)
def testNegativeIndex(self):
- "Negative index counts backwards (like python sequences)."
+ "Bodies: Negative index counts backwards (like python sequences)."
self.assert_(O.bodies[-1]==O.bodies[self.count-1])
def testErasedIterate(self):
- "Iterator over O.bodies silently skips erased bodies."
+ "Bodies: Iterator silently skips erased ones"
removed,counted=0,0
for i in range(0,10):
id=random.randint(0,self.count-1)
@@ -92,27 +92,27 @@
utils.sphere([1,1,1],.5,material=1)
])
def testShared(self):
- "shared_ptr's makes change in material immediate everywhere"
+ "Material: shared_ptr's makes change in material immediate everywhere"
O.bodies[0].mat.young=23423333
self.assert_(O.bodies[0].mat.young==O.bodies[1].mat.young)
def testSharedAfterReload(self):
- "shared_ptr's are preserved when loading from XML (using a hack in MetaBody::postProcessAttributes)"
+ "Material: shared_ptr's are preserved when saving/loading"
O.saveTmp(); O.loadTmp()
O.bodies[0].mat.young=9087438484
self.assert_(O.bodies[0].mat.young==O.bodies[1].mat.young)
def testLen(self):
- "len(O.materials)"
+ "Material: len(O.materials)"
self.assert_(len(O.materials)==2)
def testNegativeIndex(self):
- "Negative index counts backwards."
+ "Material: negative index counts backwards."
self.assert_(O.materials[-1]==O.materials[1])
def testIterate(self):
- "Iteration over O.materials"
+ "Material: iteration over O.materials"
counted=0
for m in O.materials: counted+=1
self.assert_(counted==len(O.materials))
def testAccess(self):
- "Material found by index or label; KeyError raised for invalid label."
+ "Material: find by index or label; KeyError raised for invalid label."
self.assertRaises(KeyError,lambda: O.materials['nonexistent label'])
self.assert_(O.materials['materialZero']==O.materials[0])
=== renamed file 'py/tests/cell.py' => 'py/tests/pbc.py'
--- py/tests/cell.py 2010-10-18 19:27:07 +0000
+++ py/tests/pbc.py 2010-10-19 14:57:27 +0000
@@ -14,7 +14,8 @@
from yade import *
-class TestCell(unittest.TestCase):
+class TestPBC(unittest.TestCase):
+ # prefix test names with PBC:
def setUp(self):
O.reset(); O.periodic=True;
O.cell.refSize=Vector3(2.5,2.5,3)
@@ -24,34 +25,33 @@
O.bodies.append(utils.sphere((1,1,1),.5))
self.initPos=Vector3([O.bodies[0].state.pos[i]+self.relDist[i]+self.cellDist[i]*O.cell.refSize[i] for i in (0,1,2)])
O.bodies.append(utils.sphere(self.initPos,.5))
- print O.bodies[1].state.pos
+ #print O.bodies[1].state.pos
O.bodies[1].state.vel=self.initVel
O.dt=1e-5
O.engines=[NewtonIntegrator(homotheticCellResize=2)]
O.cell.velGrad=Matrix3(0,0,0, 0,0,0, 0,0,-1)
def testHomotheticResize(self):
- "Homothetic cell resize adjusts particle velocity"
+ "PBC: homothetic cell resize adjusts particle velocity"
O.step()
s1=O.bodies[1].state
self.assertAlmostEqual(s1.vel[2],self.initVel[2]+self.initPos[2]*O.cell.velGrad[2,2])
def testScGeomIncidentVelocity(self):
+ "PBC: ScGeom computes incident velocity correctly"
+ O.dt=0 # do not change positions with dt=0 in NewtonIntegrator, but still update velocities from velGrad
O.step()
O.engines=[InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[])]
i=utils.createInteraction(0,1)
- if 0:
- print i.geom.incidentVel(i,avoidGranularRatcheting=True)
- print i.geom.incidentVel(i,avoidGranularRatcheting=False)
- print i.geom.penetrationDepth
- print i.geom.dict()
- print '=========='
- print 'v1',O.bodies[1].state.vel
- print 'shiftVel',O.cell.velGrad*O.cell.Hsize*(0,0,-10)
-
+ self.assertEqual(self.initVel,i.geom.incidentVel(i,avoidGranularRatcheting=True))
+ self.assertEqual(self.initVel,i.geom.incidentVel(i,avoidGranularRatcheting=False))
+ self.assertAlmostEqual(self.relDist[1],1-i.geom.penetrationDepth)
def testKineticEnergy(self):
+ "PBC: utils.kineticEnergy considers only fluctuation velocity, not the velocity gradient"
+ O.dt=1e-50 # with timestep almost zero, just update the velocity in the integrator
O.step() # updates velocity with homotheticCellResize
# ½(mv²+ωIω)
# #0 is still, no need to add it; #1 has zero angular velocity
# we must take self.initVel since O.bodies[1].state.vel now contains the homothetic resize which utils.kineticEnergy is supposed to compensate back
Ek=.5*O.bodies[1].state.mass*self.initVel.squaredNorm()
+ self.assertAlmostEqual(Ek,utils.kineticEnergy())
=== modified file 'py/tests/wrapper.py'
--- py/tests/wrapper.py 2010-10-18 19:27:07 +0000
+++ py/tests/wrapper.py 2010-10-19 14:57:27 +0000
@@ -22,24 +22,24 @@
def setUp(self):
pass # no setup needed for tests here
def testClassCtors(self):
- "Correct types are instantiated"
+ "Core: correct types are instantiated"
# correct instances created with Foo() syntax
for r in allClasses:
obj=eval(r)();
self.assert_(obj.__class__.__name__==r,'Failed for '+r)
def testRootDerivedCtors_attrs_few(self):
- "Class ctor's attributes"
+ "Core: class ctor's attributes"
# attributes passed when using the Foo(attr1=value1,attr2=value2) syntax
gm=Shape(wire=True); self.assert_(gm.wire==True)
def testDispatcherCtor(self):
- "Dispatcher ctors with functors"
+ "Core: dispatcher ctors with functors"
# dispatchers take list of their functors in the ctor
# same functors are collapsed in one
cld1=LawDispatcher([Law2_Dem3DofGeom_FrictPhys_CundallStrack(),Law2_Dem3DofGeom_FrictPhys_CundallStrack()]); self.assert_(len(cld1.functors)==1)
# two different make two different, right?
cld2=LawDispatcher([Law2_Dem3DofGeom_FrictPhys_CundallStrack(),Law2_Dem3DofGeom_CpmPhys_Cpm()]); self.assert_(len(cld2.functors)==2)
def testInteractionLoopCtor(self):
- "InteractionLoop special ctor"
+ "Core: InteractionLoop special ctor"
# InteractionLoop takes 3 lists
id=InteractionLoop([Ig2_Facet_Sphere_Dem3DofGeom(),Ig2_Sphere_Sphere_Dem3DofGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_Dem3DofGeom_FrictPhys_CundallStrack()],)
self.assert_(len(id.geomDispatcher.functors)==2)
@@ -47,7 +47,7 @@
self.assert_(id.physDispatcher.functors[0].__class__==Ip2_FrictMat_FrictMat_FrictPhys().__class__)
self.assert_(id.lawDispatcher.functors[0].__class__==Law2_Dem3DofGeom_FrictPhys_CundallStrack().__class__)
def testParallelEngineCtor(self):
- "ParallelEngine special ctor"
+ "Core: ParallelEngine special ctor"
pe=ParallelEngine([InsertionSortCollider(),[BoundDispatcher(),ForceResetter()]])
self.assert_(pe.slaves[0].__class__==InsertionSortCollider().__class__)
self.assert_(len(pe.slaves[1])==2)
@@ -57,9 +57,11 @@
## testing incorrect operations that should raise exceptions
##
def testWrongFunctorType(self):
+ "Core: dispatcher and functor type mismatch is detected"
# dispatchers accept only correct functors
self.assertRaises(TypeError,lambda: LawDispatcher([Bo1_Sphere_Aabb()]))
def testInvalidAttr(self):
+ 'Core: invalid attribute access raises AttributeError'
# accessing invalid attributes raises AttributeError
self.assertRaises(AttributeError,lambda: Sphere(attributeThatDoesntExist=42))
self.assertRaises(AttributeError,lambda: Sphere().attributeThatDoesntExist)
@@ -70,13 +72,13 @@
self.assertEqual(len(v1),len(v2));
for i in range(len(v1)): self.assertAlmostEqual(v1[i],v2[i],msg='Component '+str(i)+' of '+str(v1)+' and '+str(v2))
def testVector2(self):
- "Vector2 operations"
+ "Math: Vector2 operations"
v=Vector2(1,2); v2=Vector2(3,4)
self.assert_(v+v2==Vector2(4,6))
self.assert_(Vector2().UnitX.dot(Vector2().UnitY)==0)
self.assert_(Vector2().Zero.norm()==0)
def testVector3(self):
- "Vector3 operations"
+ "Math: Vector3 operations"
v=Vector3(3,4,5); v2=Vector3(3,4,5)
self.assert_(v[0]==3 and v[1]==4 and v[2]==5)
self.assert_(v.squaredNorm()==50)
@@ -88,7 +90,7 @@
self.assert_(x.dot(y)==0)
self.assert_(x.cross(y)==z)
def testQuaternion(self):
- "Quaternion operations"
+ "Math: Quaternion operations"
# construction
q1=Quaternion((0,0,1),pi/2)
q2=Quaternion(Vector3(0,0,1),pi/2)
@@ -100,7 +102,7 @@
self.assertSeqAlmostEqual(q1.toAxisAngle()[0],(0,0,1))
self.assertAlmostEqual(q1.toAxisAngle()[1],pi/2)
def testMatrix3(self):
- "Matrix3 operations"
+ "Math: Matrix3 operations"
#construction
m1=Matrix3(1,0,0,0,1,0,0,0,1)
# comparison
=== modified file 'py/wrapper/yadeWrapper.cpp'
--- py/wrapper/yadeWrapper.cpp 2010-09-30 18:00:41 +0000
+++ py/wrapper/yadeWrapper.cpp 2010-10-19 14:57:27 +0000
@@ -326,7 +326,7 @@
void dt_set(double dt){
Scene* scene=OMEGA.getScene().get();
// activate timestepper, if possible (throw exception if there is none)
- if(dt<=0){ if(!scene->timeStepperActivate(true)) /* not activated*/ throw runtime_error("No TimeStepper found in O.engines."); if(dt!=0) scene->dt=-dt; }
+ if(dt<0){ if(!scene->timeStepperActivate(true)) /* not activated*/ throw runtime_error("No TimeStepper found in O.engines."); }
else { scene->timeStepperActivate(false); scene->dt=dt; }
}
bool dynDt_get(){return OMEGA.getScene()->timeStepperActive();}
@@ -511,7 +511,7 @@
.add_property("stopAtIter",&pyOmega::stopAtIter_get,&pyOmega::stopAtIter_set,"Get/set number of iteration after which the simulation will stop.")
.add_property("time",&pyOmega::time,"Return virtual (model world) time of the simulation.")
.add_property("realtime",&pyOmega::realTime,"Return clock (human world) time the simulation has been running.")
- .add_property("dt",&pyOmega::dt_get,&pyOmega::dt_set,"Current timestep (Δt) value.\n\n* assigning zero enables dynamic Δt control via a :yref:`TimeStepper` (raises an exception if there is no :yref:`TimeStepper` among :yref:`O.engines<Omega.engines>`)\n* assigning negative value enables dynamic Δt (as in the previous case) and sets positive timestep ``O.dt=|Δt|`` (will be used until the timestepper is run and updates it)\n* assigning positive value sets Δt to that value and disables dynamic Δt (via :yref:`TimeStepper`, if there is one).\n\n:yref:`dynDt<Omega.dynDt>` can be used to query whether dynamic Δt is in use.")
+ .add_property("dt",&pyOmega::dt_get,&pyOmega::dt_set,"Current timestep (Δt) value.\n\n* assigning negative value enables dynamic Δt (by looking for a :yref:`TimeStepper` in :yref:`O.engine<Omega.engines>`) and sets positive timestep ``O.dt=|Δt|`` (will be used until the timestepper is run and updates it)\n* assigning positive value sets Δt to that value and disables dynamic Δt (via :yref:`TimeStepper`, if there is one).\n\n:yref:`dynDt<Omega.dynDt>` can be used to query whether dynamic Δt is in use.")
.add_property("dynDt",&pyOmega::dynDt_get,"Whether a :yref:`TimeStepper` is used for dynamic Δt control. See :yref:`dt<Omega.dt>` on how to enable/disable :yref:`TimeStepper`.")
.add_property("dynDtAvailable",&pyOmega::dynDtAvailable_get,"Whether a :yref:`TimeStepper` is amongst :yref:`O.engines<Omega.engines>`, activated or not.")
.def("load",&pyOmega::load,"Load simulation from file.")
