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 _______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
