yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01451
[svn] r1853 - in trunk: pkg/common pkg/common/Engine/StandAloneEngine pkg/dem/Engine/DeusExMachina scripts/test
Author: eudoxos
Date: 2009-07-10 21:16:01 +0200 (Fri, 10 Jul 2009)
New Revision: 1853
Modified:
trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp
trunk/pkg/common/SConscript
trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp
trunk/scripts/test/collider-stride-triax.py
Log:
1. Make NewtonsDampedLaw compute maxVlocitySq at every step anew.
2. Make InsertionSortCollider find NewtonsDampedLaw and get max speed from there. Together with InsertionSortColider::sweepLength, this makes stride adapted automatically to maxVelocitySq; this will described at http://yade.wikia.com/wiki/Insertion_Sort_Collider_Stride.
3. Update test script for that.
Modified: trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp 2009-07-10 09:52:33 UTC (rev 1852)
+++ trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp 2009-07-10 19:16:01 UTC (rev 1853)
@@ -5,6 +5,7 @@
#include<yade/core/Interaction.hpp>
#include<yade/core/InteractionContainer.hpp>
#include<yade/pkg-common/BoundingVolumeMetaEngine.hpp>
+#include<yade/pkg-dem/NewtonsDampedLaw.hpp>
#include<algorithm>
#include<vector>
@@ -65,16 +66,31 @@
// activated if number of bodies changes (hence need to refresh collision information)
// or the time of scheduled run already came, or we were never scheduled yet
if(stride<=1) return true;
- bool ret=XX.size()!=2*rb->bodies->size() || scheduledRun<0 || rb->simulationTime>=scheduledRun;
- // we wouldn't run in this step; just delete pending interactions
- // this could be done in ::action, but it would make the call counters not reflect the stride
+ if(sweepLength<=0) return true;
+ if(!newton) return true; // we wouldn't be able to find the max velocity
+ bool ret=rb->simulationTime>=scheduledRun ||
+ // if the max velocity is bigger than the one that we used for bb computation last time
+ // and the distance bodies would travel with this bigger velocity since last run (rb->simulationTime-lastRun)
+ // would be the same or greater than the one that would be traveled with the original velocity
+ // over the stride time (scheduledRun-lastRun)
+ (newton->maxVelocitySq<0) || // no valid data about max velocity, run always
+ (sweepVelocity==0 && newton->maxVelocitySq>0) ||
+ (sweepVelocity<sqrt(newton->maxVelocitySq) && rb->simulationTime-lastRun>=(sweepVelocity/sqrt(newton->maxVelocitySq)/* we know maxVelocitySq>0 from the first condition */)*(scheduledRun-lastRun)) ||
+ // number of bodies changed
+ XX.size()!=2*rb->bodies->size() ||
+ // we've never run yet (this should never happen as per if(!newton) above)
+ scheduledRun<0;
+ // we wouldn't run in this step; in that case, just delete pending interactions
+ // this is done in ::action normally, but it would make the call counters not reflect the stride
if(!ret) rb->interactions->erasePending(*this);
return ret;
}
#endif
void InsertionSortCollider::action(MetaBody* rb){
- // timingDeltas->start();
+ #ifdef ISC_TIMING
+ timingDeltas->start();
+ #endif
size_t nBodies=rb->bodies->size();
InteractionContainer* interactions=rb->interactions.get();
@@ -107,28 +123,44 @@
// get the BoundingVolumeMetaEngine and turn it off; we will call it ourselves
if(!boundDispatcher){
FOREACH(shared_ptr<Engine>& e, rb->engines){ boundDispatcher=dynamic_pointer_cast<BoundingVolumeMetaEngine>(e); if(boundDispatcher) break; }
- if(!boundDispatcher){ LOG_FATAL("Unable to locate BoundingVolumeMetaEngine within engines, aborting."); throw runtime_error("Explanation above"); }
+ if(!boundDispatcher){ LOG_FATAL("Unable to locate BoundingVolumeMetaEngine within engines, aborting."); abort(); }
boundDispatcher->activated=false; // deactive the engine, we will call it ourselves from now (just when needed)
}
+ if(sweepLength>0){
+ // get NewtonsDampedLaw, to ask for the maximum velocity value
+ if(!newton){
+ FOREACH(shared_ptr<Engine>& e, rb->engines){ newton=dynamic_pointer_cast<NewtonsDampedLaw>(e); if(newton) break; }
+ if(!newton){ LOG_FATAL("Unable to locate NewtonsDampedLaw within engines, aborting."); abort(); }
+ }
+ }
#endif
- // timingDeltas->checkpoint("init");
+ ISC_CHECKPOINT("init");
#ifdef COLLIDE_STRIDED
- // FIXME: should be able to adapt stride based on the potential_interaction_count_increase/stride_speedup tradeoff
- // this depends on the packing density, for instance, maximum velocities/accels etc.
- if(stride>1 && sweepTimeFactor<1 && sweepVelocity<=0){ LOG_WARN("Stride is "<<stride<<", but no sweeping effective!! Setting stride back to 1."); stride=1; }
- if(stride>1){
- //schedule next run
- scheduledRun=rb->simulationTime+rb->dt*(stride-.5); // -.5 to avoid rounding issues
- if(sweepTimeFactor>=1) boundDispatcher->sweepTime=rb->dt*stride*sweepTimeFactor;
- if(sweepVelocity>0) boundDispatcher->sweepDist=rb->dt*stride*sweepVelocity;
+ if(sweepLength>0){
+ if(newton->maxVelocitySq>=0){ // non-negative, i.e. a really computed value
+ // compute new stride value
+ assert(sweepFactor>1.);
+ sweepVelocity=sqrt(newton->maxVelocitySq)*sweepFactor;
+ if(sweepVelocity>0) {
+ stride=max(1,int((sweepLength/sweepVelocity)/rb->dt));
+ boundDispatcher->sweepDist=rb->dt*(stride-1)*sweepVelocity;
+ } else { // no motion
+ stride=1000; // shouldn't this be some saner value? Infinity? How to decide?
+ boundDispatcher->sweepDist=0; // nothing moves, no need to make bboxes larger
+ }
+ scheduledRun=rb->simulationTime+rb->dt*(stride-.5); // -.5 to avoid rounding issues
+ } else { /* no valid data yet, run next time again */ boundDispatcher->sweepDist=0; stride=1; scheduledRun=rb->simulationTime+rb->dt; }
+ LOG_DEBUG(rb->simulationTime<<"s: stride adapted to "<<stride<<"; sweepVelocity="<<sweepVelocity<<", maxVelocity="<<sqrt(newton->maxVelocitySq)<<", sweepDist="<<boundDispatcher->sweepDist);
+ newton->maxVelocitySq=-1; // reset to invalid value again
} else { scheduledRun=-1; boundDispatcher->sweepTime=-1; boundDispatcher->sweepDist=0; }
boundDispatcher->action(rb);
#endif
- // timingDeltas->checkpoint("bound");
+ ISC_CHECKPOINT("bound");
+
// copy bounds along given axis into our arrays
for(size_t i=0; i<2*nBodies; i++){
const body_id_t& idXX=XX[i].id; const body_id_t& idYY=YY[i].id; const body_id_t& idZZ=ZZ[i].id;
@@ -144,11 +176,11 @@
else{ const Vector3r& pos=Body::byId(idXX,rb)->physicalParameters->se3.position; memcpy(&minima[3*idXX],pos,3*sizeof(Real)); memcpy(&maxima[3*idXX],pos,3*sizeof(Real)); }
}
}
+ ISC_CHECKPOINT("copy");
- // timingDeltas->checkpoint("copy");
-
// process interactions that the constitutive law asked to be erased
- interactions->erasePending(*this);
+ interactions->erasePending(*this);
+ ISC_CHECKPOINT("erase");
// sort
if(!doInitSort && !sortThenCollide){
@@ -159,13 +191,9 @@
if(doInitSort){
// the initial sort is in independent in 3 dimensions, may be run in parallel
// it seems that there is no time gain running this in parallel, though
- #pragma omp parallel sections
{
- #pragma omp section
std::sort(XX.begin(),XX.end());
- #pragma omp section
std::sort(YY.begin(),YY.end());
- #pragma omp section
std::sort(ZZ.begin(),ZZ.end());
}
} else { // sortThenCollide
@@ -205,5 +233,5 @@
}
}
}
- // timingDeltas->checkpoint("sort&collide");
+ ISC_CHECKPOINT("sort&collide");
}
Modified: trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp 2009-07-10 09:52:33 UTC (rev 1852)
+++ trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp 2009-07-10 19:16:01 UTC (rev 1853)
@@ -15,7 +15,18 @@
*/
#define COLLIDE_STRIDED
+
+// #define this macro to enable timing whichin this engine
+//#define ISC_TIMING
+
+#ifdef ISC_TIMING
+ #define ISC_CHECKPOINT(cpt) timingDeltas->checkPoint(cpt)
+#else
+ #define ISC_CHECKPOINT(cpt)
+#endif
+
class BoundingVolumeMetaEngine;
+class NewtonsDampedLaw;
class InsertionSortCollider: public Collider{
//! struct for storing bounds of bodies
@@ -33,16 +44,25 @@
#ifdef COLLIDE_STRIDED
// keep this dispatcher and call it ourselves as needed
shared_ptr<BoundingVolumeMetaEngine> boundDispatcher;
- // interval at which we will run; if <=1, we run always (as usual). 0 by default.
+ // we need this to find out about current maxVelocitySq
+ shared_ptr<NewtonsDampedLaw> newton;
+ // interval at which we will run; if 0, we run always (as usual).
int stride;
- // virtual time when we have to run the next time
- Real scheduledRun;
- //! If >1 and using stride, sweep time will be multiplied by this number; it should be >=1. to accomodate non-linearities in the system.
- /// If deactivated (-1 by default), current-velocity-based sweeping will not be enabled.
- Real sweepTimeFactor;
- //! If >0 and using stride, all bodies will be swept as if having this velocity.
- /// It should be the maximum (predicted or measured) velocity increased by some safety margin, otherwise bodies may get out of their AABB. Deactivated (-1) by default.
+ //! virtual time when we were run last time
+ Real lastRun;
+ //! virtual time when we have to run the next time
+ Real scheduledRun;
+ /// Absolute length that will be added to bounding boxes at each side; it should be something like 1/5 of typical grain radius
+ /// this value is used to adapt stride; if too large, stride will be big, but the ratio of potential-only interactions will be very big,
+ /// thus slowing down collider & interaction loops significantly (remember: O(addLength^3))
+ /// If non-positive, collider runs always, without stride adaptivity
+ Real sweepLength;
+ //! Current sweeping velocity; computed from maxVelocity*sweepVelocityFactor
Real sweepVelocity;
+ //! Overestimation factor for the sweep velocity; must be >=1.0.
+ /// Has no influence on sweepLength, only on the computed stride.
+ /// Default 1.05
+ Real sweepFactor;
#endif
//! storage for bounds
std::vector<Bound> XX,YY,ZZ;
@@ -72,14 +92,18 @@
InsertionSortCollider():
#ifdef COLLIDE_STRIDED
- stride(0), scheduledRun(-1), sweepTimeFactor(-1), sweepVelocity(-1),
+ stride(0), lastRun(-1), scheduledRun(-1), sweepLength(-1), sweepVelocity(-1), sweepFactor(1.05),
#endif
- sortAxis(0), sortThenCollide(false){ /* timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas); */ }
+ sortAxis(0), sortThenCollide(false){
+ #ifdef ISC_TIMING
+ timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
+ #endif
+ }
virtual void action(MetaBody*);
REGISTER_CLASS_AND_BASE(InsertionSortCollider,Collider);
REGISTER_ATTRIBUTES(Collider,(sortAxis)(sortThenCollide)
#ifdef COLLIDE_STRIDED
- (stride)(scheduledRun)(sweepTimeFactor)(sweepVelocity)
+ (stride)(lastRun)(scheduledRun)(sweepLength)(sweepFactor)(sweepVelocity)
#endif
);
DECLARE_LOGGER;
Modified: trunk/pkg/common/SConscript
===================================================================
--- trunk/pkg/common/SConscript 2009-07-10 09:52:33 UTC (rev 1852)
+++ trunk/pkg/common/SConscript 2009-07-10 19:16:01 UTC (rev 1853)
@@ -143,7 +143,7 @@
env.SharedLibrary('SpheresFactory',['Engine/StandAloneEngine/SpheresFactory.cpp'],
LIBS=env['LIBS']+['AABB','InteractingSphere','Facet','Sphere','BodyMacroParameters','InteractionGeometryMetaEngine']),
env.SharedLibrary('SpatialQuickSortCollider',['Engine/StandAloneEngine/SpatialQuickSortCollider.cpp']),
- env.SharedLibrary('InsertionSortCollider',['Engine/StandAloneEngine/InsertionSortCollider.cpp'],LIBS=env['LIBS']+['BoundingVolumeMetaEngine']),
+ env.SharedLibrary('InsertionSortCollider',['Engine/StandAloneEngine/InsertionSortCollider.cpp'],LIBS=env['LIBS']+['BoundingVolumeMetaEngine','NewtonsDampedLaw']),
env.SharedLibrary('PersistentSAPCollider',['Engine/StandAloneEngine/PersistentSAPCollider.cpp']),
env.SharedLibrary('DistantPersistentSAPCollider',['Engine/StandAloneEngine/DistantPersistentSAPCollider.cpp']),
env.SharedLibrary('PhysicalActionContainerReseter',['Engine/StandAloneEngine/PhysicalActionContainerReseter.cpp']),
Modified: trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp 2009-07-10 09:52:33 UTC (rev 1852)
+++ trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp 2009-07-10 19:16:01 UTC (rev 1853)
@@ -25,6 +25,7 @@
{
ncb->bex.sync();
Real dt=Omega::instance().getTimeStep();
+ maxVelocitySq=-1;
FOREACH(const shared_ptr<Body>& b, *ncb->bodies){
// clump members are non-dynamic; they skip the rest of loop once their forces are properly taken into account, however
Modified: trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp 2009-07-10 09:52:33 UTC (rev 1852)
+++ trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp 2009-07-10 19:16:01 UTC (rev 1853)
@@ -38,10 +38,10 @@
public:
///damping coefficient for Cundall's non viscous damping
Real damping;
- /// store square of max. velocity, for informative purposes
+ /// store square of max. velocity, for informative purposes; computed again at every step
Real maxVelocitySq;
virtual void applyCondition(MetaBody *);
- NewtonsDampedLaw(): damping(0.2), maxVelocitySq(0){}
+ NewtonsDampedLaw(): damping(0.2), maxVelocitySq(-1){}
REGISTER_ATTRIBUTES(DeusExMachina,(damping)(maxVelocitySq));
REGISTER_CLASS_AND_BASE(NewtonsDampedLaw,DeusExMachina);
};
Modified: trunk/scripts/test/collider-stride-triax.py
===================================================================
--- trunk/scripts/test/collider-stride-triax.py 2009-07-10 09:52:33 UTC (rev 1852)
+++ trunk/scripts/test/collider-stride-triax.py 2009-07-10 19:16:01 UTC (rev 1853)
@@ -4,30 +4,19 @@
import os.path
loadFrom='/tmp/triax.xml'
if not os.path.exists(loadFrom):
- TriaxialTest(numberOfGrains=8000,fast=fast,noFiles=True).generate(loadFrom)
+ TriaxialTest(numberOfGrains=8000,fast=True,noFiles=True).generate(loadFrom)
O.load(loadFrom)
log.setLevel('TriaxialCompressionEngine',log.WARN) # shut up
+log.setLevel('InsertionSortCollider',log.DEBUG) # shut up
collider=utils.typedEngine('InsertionSortCollider')
newton=utils.typedEngine('NewtonsDampedLaw')
# use striding; say "if 0:" to disable striding and compare to regular runs
if 1:
- collider['stride']=4
- collider['sweepVelocity']=16
- # try to update maximum velocity dynamically:
- if 1:
- O.engines=O.engines+[PeriodicPythonRunner(command='adjustMaxVelocity()',iterPeriod=50)]
+ # length by which bboxes will be made larger
+ collider['sweepLength']=.05*O.bodies[100].shape['radius']
-def adjustMaxVelocity():
- newMax=sqrt(newton['maxVelocitySq']); newton['maxVelocitySq']=0
- collider['sweepVelocity']=newMax
- # try to adjust stride here; the 4 and 16 are original (reference) values for stride and sweepVelocity
- # if this were inverse-proportional (without sqrt), stride is getting too high
- collider['stride']=4*int(sqrt(16/newMax))
- print 'step %d; new max velocity %g, stride %d'%(O.iter,newMax,collider['stride'])
-
-
O.step() # filter out initialization
O.timingEnabled=True
totalTime=0
Follow ups
