Author: eudoxos Date: 2009-03-29 17:17:31 +0200 (Sun, 29 Mar 2009) New Revision: 1736
Added: trunk/examples/triax-perf/ trunk/examples/triax-perf/triax-perf.py trunk/examples/triax-perf/triax-perf.table Modified: trunk/core/yade.cpp trunk/extra/Brefcom.cpp trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.hpp trunk/pkg/dem/PreProcessor/TriaxialTest.cpp trunk/pkg/dem/PreProcessor/TriaxialTest.hpp trunk/pkg/dem/SConscript trunk/pkg/fem/Engine/StandAloneEngine/FEMLaw.cpp trunk/pkg/mass-spring/Engine/StandAloneEngine/MassSpringLaw.cpp Log: 1. Create ef2_Spheres_Elastic_ElasticLaw that has the ElasticContactLaw algorithm 2. ElasticContactLaw now merely calls ef2_Spheres_Elastic_ElasticLaw. 3. TriaxialTest::parallel controls whether to use InteractionDispatchers or not. 4. Added examples/triax-perf to show the impact of that. Also see http://yade.wikia.com/wiki/Triaxial_Test_Parallel 5. Fix a few compilation issues. 6. Flush stdout before exiting from main, so that in case of crash upon exit (log4cxx?) we have all messages output. Modified: trunk/core/yade.cpp =================================================================== --- trunk/core/yade.cpp 2009-03-29 10:34:26 UTC (rev 1735) +++ trunk/core/yade.cpp 2009-03-29 15:17:31 UTC (rev 1736) @@ -303,6 +303,7 @@ #endif LOG_INFO("Yade: normal exit."); + fflush(stdout); // in case of crash at exit, logs will not disappear return ok; } Added: trunk/examples/triax-perf/triax-perf.py =================================================================== --- trunk/examples/triax-perf/triax-perf.py 2009-03-29 10:34:26 UTC (rev 1735) +++ trunk/examples/triax-perf/triax-perf.py 2009-03-29 15:17:31 UTC (rev 1736) @@ -0,0 +1,23 @@ +# Performance test for running +# +# 1. Regular TriaxialTest with 3 independent dispatchers (geom, phys, constitutive law) +# 2. TriaxialTest with InteractionDispatchers (common loop and functor cache) +# +# Run the test like this: +# +# yade-trunk-opt-multi -j1 triax-perf.table triax-perf.py +# +# The -j1 ensures that only 1 job will run at time +# (even if other cores are free, access to memory is limiting if running multiple jobs at time) +# +# You have to collect the results by hand from log files. +# +utils.readParamsFromTable(parallel=False,noTableOk=True) +p=Preprocessor('TriaxialTest',{'numberOfGrains':8000,'parallel':parallel}).load() +O.run(10,True) # filter out initialization +O.timingEnabled=True +O.run(1000,True) +from yade import timing +timing.stats() +print 'BexContainer synced %d times'%(O.bexSyncCount) + Added: trunk/examples/triax-perf/triax-perf.table =================================================================== --- trunk/examples/triax-perf/triax-perf.table 2009-03-29 10:34:26 UTC (rev 1735) +++ trunk/examples/triax-perf/triax-perf.table 2009-03-29 15:17:31 UTC (rev 1736) @@ -0,0 +1,11 @@ +!OMP_NUM_THREADS parallel description +1 False ser1 +2 False ser2 +3 False ser3 +4 False ser4 +5 False ser5 +1 True par1 +2 True par2 +3 True par3 +4 True par4 +5 True par5 Modified: trunk/extra/Brefcom.cpp =================================================================== --- trunk/extra/Brefcom.cpp 2009-03-29 10:34:26 UTC (rev 1735) +++ trunk/extra/Brefcom.cpp 2009-03-29 15:17:31 UTC (rev 1736) @@ -233,8 +233,7 @@ } -void BrefcomLaw::action(MetaBody* _rootBody){ - rootBody=_rootBody; +void BrefcomLaw::action(MetaBody* rootBody){ if(!functor) functor=shared_ptr<ef2_Spheres_Brefcom_BrefcomLaw>(new ef2_Spheres_Brefcom_BrefcomLaw); functor->logStrain=logStrain; FOREACH(const shared_ptr<Interaction>& I, *rootBody->interactions){ Modified: trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp =================================================================== --- trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp 2009-03-29 10:34:26 UTC (rev 1735) +++ trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp 2009-03-29 15:17:31 UTC (rev 1736) @@ -20,31 +20,36 @@ FOREACH(shared_ptr<Interaction> I, *rootBody->interactions){ #endif #ifdef DISPATCH_CACHE - shared_ptr<Body>& b1=(*rootBody->bodies)[I->getId1()]; - shared_ptr<Body>& b2=(*rootBody->bodies)[I->getId2()]; + const shared_ptr<Body>& b1_=Body::byId(I->getId1(),rootBody); + const shared_ptr<Body>& b2_=Body::byId(I->getId2(),rootBody); + // we know there is no geometry functor already, take the short path if(!I->functorCache.geomExists) { I->isReal=false; continue; } // no interaction geometry for either of bodies; no interaction possible - if(!b1->interactingGeometry || !b2->interactingGeometry) { I->isReal=false; continue; } + if(!b1_->interactingGeometry || !b2_->interactingGeometry) { I->isReal=false; continue; } + bool swap=false; // InteractionGeometryMetaEngine if(!I->functorCache.geom || !I->functorCache.phys){ - bool swap=false; - I->functorCache.geom=geomDispatcher->getFunctor2D(b1->interactingGeometry,b2->interactingGeometry,swap); + I->functorCache.geom=geomDispatcher->getFunctor2D(b1_->interactingGeometry,b2_->interactingGeometry,swap); // returns NULL ptr if no functor exists; remember that and shortcut if(!I->functorCache.geom) { I->functorCache.geomExists=false; continue; } - // arguments for the geom functor are in the reverse order (dispatcher would normally call goReverse). - // we don't remember the fact that is reverse, so we swap bodies within the interaction - // and can call go in all cases - if(swap) { I->swapOrder(); b1=(*rootBody->bodies)[I->getId1()]; b2=(*rootBody->bodies)[I->getId2()]; } } + // arguments for the geom functor are in the reverse order (dispatcher would normally call goReverse). + // we don't remember the fact that is reverse, so we swap bodies within the interaction + // and can call go in all cases + if(swap){I->swapOrder();} + // body pointers must be updated, in case we swapped + const shared_ptr<Body>& b1=Body::byId(I->getId1()); + const shared_ptr<Body>& b2=Body::byId(I->getId2()); + assert(I->functorCache.geom); I->isReal=I->functorCache.geom->go(b1->interactingGeometry,b2->interactingGeometry,b1->physicalParameters->se3, b2->physicalParameters->se3,I); if(!I->isReal) continue; // InteractionPhysicsMetaEngine if(!I->functorCache.phys){ - bool swap=false; I->functorCache.phys=physDispatcher->getFunctor2D(b1->physicalParameters,b2->physicalParameters,swap); + I->functorCache.phys=physDispatcher->getFunctor2D(b1->physicalParameters,b2->physicalParameters,swap); assert(!swap); // InteractionPhysicsEngineUnits are symmetric } assert(I->functorCache.phys); @@ -54,7 +59,7 @@ // populating constLaw cache must be done after geom and physics dispatchers have been called, since otherwise the interaction // would not have interactionGeometry and interactionPhysics yet. if(!I->functorCache.constLaw){ - bool swap=false; I->functorCache.constLaw=constLawDispatcher->getFunctor2D(I->interactionGeometry,I->interactionPhysics,swap); + I->functorCache.constLaw=constLawDispatcher->getFunctor2D(I->interactionGeometry,I->interactionPhysics,swap); assert(!swap); // reverse call would make no sense, as the arguments are of different types } assert(I->functorCache.constLaw); Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp =================================================================== --- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-03-29 10:34:26 UTC (rev 1735) +++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-03-29 15:17:31 UTC (rev 1736) @@ -19,6 +19,8 @@ #include<yade/extra/Shop.hpp> +YADE_PLUGIN("ElasticContactLaw2","ef2_Spheres_Elastic_ElasticLaw","ElasticContactLaw"); + ElasticContactLaw2::ElasticContactLaw2(){ Shop::Bex::initCache(); isCohesive=true; @@ -58,12 +60,18 @@ -ElasticContactLaw::ElasticContactLaw() : InteractionSolver() , actionForce(new Force) , actionMomentum(new Momentum) + +ElasticContactLaw::ElasticContactLaw() : InteractionSolver() +#ifndef BEX_CONTAINER + , actionForce(new Force) , actionMomentum(new Momentum) +#endif { sdecGroupMask=1; momentRotationLaw = true; - actionForceIndex = actionForce->getClassIndex(); - actionMomentumIndex = actionMomentum->getClassIndex(); + #ifndef BEX_CONTAINER + actionForceIndex = actionForce->getClassIndex(); + actionMomentumIndex = actionMomentum->getClassIndex(); + #endif #ifdef SCG_SHEAR useShear=false; #endif @@ -81,34 +89,37 @@ } -void ElasticContactLaw::action(MetaBody* ncb) +void ElasticContactLaw::action(MetaBody* rootBody) { - shared_ptr<BodyContainer>& bodies = ncb->bodies; + if(!functor) functor=shared_ptr<ef2_Spheres_Elastic_ElasticLaw>(new ef2_Spheres_Elastic_ElasticLaw); + functor->momentRotationLaw=momentRotationLaw; + functor->sdecGroupMask=sdecGroupMask; + #ifdef SCG_SHEAR + functor->useShear=useShear; + #endif + FOREACH(const shared_ptr<Interaction>& I, *rootBody->interactions){ + if(!I->isReal) continue; + #ifdef YADE_DEBUG + // these checks would be redundant in the functor (ConstitutiveLawDispatcher does that already) + if(!dynamic_cast<SpheresContactGeometry*>(I->interactionGeometry.get()) || !dynamic_cast<ElasticContactInteraction*>(I->interactionPhysics.get())) continue; + #endif + functor->go(I->interactionGeometry, I->interactionPhysics, I.get(), rootBody); + } +} +void ef2_Spheres_Elastic_ElasticLaw::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, MetaBody* ncb){ Real dt = Omega::instance().getTimeStep(); -/// Non Permanents Links /// - - InteractionContainer::iterator ii = ncb->transientInteractions->begin(); - InteractionContainer::iterator iiEnd = ncb->transientInteractions->end(); - for( ; ii!=iiEnd ; ++ii ) - { - if ((*ii)->isReal) - { - const shared_ptr<Interaction>& contact = *ii; - int id1 = contact->getId1(); - int id2 = contact->getId2(); - - if( !( (*bodies)[id1]->getGroupMask() & (*bodies)[id2]->getGroupMask() & sdecGroupMask) ) continue; - - SpheresContactGeometry* currentContactGeometry= YADE_CAST<SpheresContactGeometry*>(contact->interactionGeometry.get()); - ElasticContactInteraction* currentContactPhysics = YADE_CAST<ElasticContactInteraction*> (contact->interactionPhysics.get()); - if((!currentContactGeometry)||(!currentContactPhysics)) continue; + int id1 = contact->getId1(), id2 = contact->getId2(); + // FIXME: mask handling should move to ConstitutiveLaw itself, outside the functors + if( !(Body::byId(id1,ncb)->getGroupMask() & Body::byId(id2,ncb)->getGroupMask() & sdecGroupMask) ) return; + SpheresContactGeometry* currentContactGeometry= static_cast<SpheresContactGeometry*>(ig.get()); + ElasticContactInteraction* currentContactPhysics = static_cast<ElasticContactInteraction*>(ip.get()); // delete interaction where spheres don't touch - if(currentContactGeometry->penetrationDepth<0){ (*ii)->isReal=false; continue; } + if(currentContactGeometry->penetrationDepth<0){ contact->isReal=false; return; } - BodyMacroParameters* de1 = YADE_CAST<BodyMacroParameters*>((*bodies)[id1]->physicalParameters.get()); - BodyMacroParameters* de2 = YADE_CAST<BodyMacroParameters*>((*bodies)[id2]->physicalParameters.get()); + BodyMacroParameters* de1 = YADE_CAST<BodyMacroParameters*>(Body::byId(id1,ncb)->physicalParameters.get()); + BodyMacroParameters* de2 = YADE_CAST<BodyMacroParameters*>(Body::byId(id1,ncb)->physicalParameters.get()); Vector3r& shearForce = currentContactPhysics->shearForce; @@ -117,7 +128,6 @@ Real un=currentContactGeometry->penetrationDepth; currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal; - // the same as under #else, but refactored #ifdef SCG_SHEAR if(useShear){ currentContactGeometry->updateShear(de1,de2,dt); @@ -156,11 +166,6 @@ static_cast<Momentum*>( ncb->physicalActions->find( id1 , actionMomentumIndex ).get() )->momentum -= _c1x.Cross(f); static_cast<Momentum*>( ncb->physicalActions->find( id2 , actionMomentumIndex ).get() )->momentum += _c2x.Cross(f); #endif - currentContactPhysics->prevNormal = currentContactGeometry->normal; - } - } } - -YADE_PLUGIN(); Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.hpp =================================================================== --- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.hpp 2009-03-29 10:34:26 UTC (rev 1735) +++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.hpp 2009-03-29 15:17:31 UTC (rev 1736) @@ -13,6 +13,7 @@ // only to see whether SCG_SHEAR is defined, may be removed in the future #include<yade/pkg-dem/SpheresContactGeometry.hpp> +#include<yade/pkg-common/ConstitutiveLaw.hpp> #include <set> #include <boost/tuple/tuple.hpp> @@ -36,14 +37,40 @@ }; REGISTER_SERIALIZABLE(ElasticContactLaw2); +class ef2_Spheres_Elastic_ElasticLaw: public ConstitutiveLaw{ + public: + virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody); + int sdecGroupMask; + bool momentRotationLaw; + #ifdef SCG_SHEAR + bool useShear; + #endif + ef2_Spheres_Elastic_ElasticLaw(): sdecGroupMask(1), momentRotationLaw(true) + #ifdef SCG_SHEAR + , useShear(false) + #endif + {} + FUNCTOR2D(SpheresContactGeometry,ElasticContactInteraction); + REGISTER_CLASS_AND_BASE(ef2_Spheres_Elastic_ElasticLaw,ConstitutiveLaw); + REGISTER_ATTRIBUTES(ConstitutiveLaw,(sdecGroupMask)(momentRotationLaw) + #ifdef SCG_SHEAR + (useShear) + #endif + ); +}; +REGISTER_SERIALIZABLE(ef2_Spheres_Elastic_ElasticLaw); + class ElasticContactLaw : public InteractionSolver { /// Attributes private : + #ifndef BEX_CONTAINER shared_ptr<PhysicalAction> actionForce; shared_ptr<PhysicalAction> actionMomentum; int actionForceIndex; int actionMomentumIndex; + NEEDS_BEX("Force","Momentum"); + #endif public : int sdecGroupMask; @@ -55,10 +82,11 @@ ElasticContactLaw(); void action(MetaBody*); + shared_ptr<ef2_Spheres_Elastic_ElasticLaw> functor; + protected : void registerAttributes(); - NEEDS_BEX("Force","Momentum"); REGISTER_CLASS_NAME(ElasticContactLaw); REGISTER_BASE_CLASS_NAME(InteractionSolver); }; Modified: trunk/pkg/dem/PreProcessor/TriaxialTest.cpp =================================================================== --- trunk/pkg/dem/PreProcessor/TriaxialTest.cpp 2009-03-29 10:34:26 UTC (rev 1735) +++ trunk/pkg/dem/PreProcessor/TriaxialTest.cpp 2009-03-29 15:17:31 UTC (rev 1736) @@ -64,8 +64,11 @@ #include<yade/pkg-common/InteractionVecSet.hpp> #include<yade/pkg-common/PhysicalActionVectorVector.hpp> +#include<yade/pkg-common/InteractionDispatchers.hpp> + #include<yade/extra/Shop.hpp> + #include <boost/filesystem/convenience.hpp> #include <boost/lexical_cast.hpp> #include <boost/numeric/conversion/bounds.hpp> @@ -162,6 +165,10 @@ isotropicCompaction=false; fixedPorosity = 1; + + #ifdef BEX_CONTAINER + parallel=false; + #endif @@ -232,8 +239,9 @@ REGISTER_ATTRIBUTE(isotropicCompaction); REGISTER_ATTRIBUTE(fixedPorosity); REGISTER_ATTRIBUTE(fixedBoxDims); - - + #ifdef BEX_CONTAINER + REGISTER_ATTRIBUTE(parallel); + #endif } @@ -253,7 +261,7 @@ //rootBody->transientInteractions = shared_ptr<InteractionContainer>(new InteractionHashMap); - rootBody->bodies = shared_ptr<BodyContainer>(new BodyRedirectionVector); + //rootBody->bodies = shared_ptr<BodyContainer>(new BodyRedirectionVector); createActors(rootBody); @@ -637,9 +645,34 @@ // rootBody->engines.push_back(sdecTimeStepper); rootBody->engines.push_back(boundingVolumeDispatcher); rootBody->engines.push_back(shared_ptr<Engine>(new PersistentSAPCollider)); - rootBody->engines.push_back(interactionGeometryDispatcher); - rootBody->engines.push_back(interactionPhysicsDispatcher); - rootBody->engines.push_back(elasticContactLaw); + #ifdef BEX_CONTAINER + if(parallel){ + #if 1 + shared_ptr<InteractionDispatchers> ids(new InteractionDispatchers); + ids->geomDispatcher=interactionGeometryDispatcher; + ids->physDispatcher=interactionPhysicsDispatcher; + ids->constLawDispatcher=shared_ptr<ConstitutiveLawDispatcher>(new ConstitutiveLawDispatcher); + shared_ptr<ef2_Spheres_Elastic_ElasticLaw> see(new ef2_Spheres_Elastic_ElasticLaw); + see->sdecGroupMask=elasticContactLaw->sdecGroupMask; + ids->constLawDispatcher->add(see); + rootBody->engines.push_back(ids); + #else + rootBody->engines.push_back(interactionGeometryDispatcher); + rootBody->engines.push_back(interactionPhysicsDispatcher); + shared_ptr<ConstitutiveLawDispatcher> cld(new ConstitutiveLawDispatcher); + shared_ptr<ef2_Spheres_Elastic_ElasticLaw> see(new ef2_Spheres_Elastic_ElasticLaw); + see->sdecGroupMask=elasticContactLaw->sdecGroupMask; + cld->add(see); + rootBody->engines.push_back(cld); + #endif + } else { + #endif + rootBody->engines.push_back(interactionGeometryDispatcher); + rootBody->engines.push_back(interactionPhysicsDispatcher); + rootBody->engines.push_back(elasticContactLaw); + #ifdef BEX_CONTAINER + } + #endif //rootBody->engines.push_back(stiffnesscounter); //rootBody->engines.push_back(stiffnessMatrixTimeStepper); Modified: trunk/pkg/dem/PreProcessor/TriaxialTest.hpp =================================================================== --- trunk/pkg/dem/PreProcessor/TriaxialTest.hpp 2009-03-29 10:34:26 UTC (rev 1735) +++ trunk/pkg/dem/PreProcessor/TriaxialTest.hpp 2009-03-29 15:17:31 UTC (rev 1736) @@ -106,6 +106,11 @@ //!flag to choose an isotropic compaction until a fixed porosity choosing a same translation speed for the six walls ,isotropicCompaction; + #ifdef BEX_CONTAINER + //! Generate parallel simulation, if it is supported + bool parallel; + #endif + int recordIntervalIter ,timeStepUpdateInterval ,timeStepOutputInterval Modified: trunk/pkg/dem/SConscript =================================================================== --- trunk/pkg/dem/SConscript 2009-03-29 10:34:26 UTC (rev 1735) +++ trunk/pkg/dem/SConscript 2009-03-29 15:17:31 UTC (rev 1736) @@ -736,6 +736,8 @@ 'TriaxialStateRecorder', 'PositionOrientationRecorder', 'Shop', + 'ConstitutiveLawDispatcher', + 'InteractionDispatchers', 'NewtonsDampedLaw']), env.SharedLibrary('CohesiveTriaxialTest', Modified: trunk/pkg/fem/Engine/StandAloneEngine/FEMLaw.cpp =================================================================== --- trunk/pkg/fem/Engine/StandAloneEngine/FEMLaw.cpp 2009-03-29 10:34:26 UTC (rev 1735) +++ trunk/pkg/fem/Engine/StandAloneEngine/FEMLaw.cpp 2009-03-29 15:17:31 UTC (rev 1736) @@ -78,7 +78,7 @@ #ifdef BEX_CONTAINER fem->bex.addForce(femTet->ids[i],force); #else - static_cast<Force*>( physicalActions->find( femTet->ids[i] , actionForce ->getClassIndex() ).get() )->force += force; + static_cast<Force*>(fem->physicalActions->find( femTet->ids[i] , actionForce ->getClassIndex() ).get() )->force += force; #endif } Modified: trunk/pkg/mass-spring/Engine/StandAloneEngine/MassSpringLaw.cpp =================================================================== --- trunk/pkg/mass-spring/Engine/StandAloneEngine/MassSpringLaw.cpp 2009-03-29 10:34:26 UTC (rev 1735) +++ trunk/pkg/mass-spring/Engine/StandAloneEngine/MassSpringLaw.cpp 2009-03-29 15:17:31 UTC (rev 1736) @@ -72,8 +72,8 @@ massSpring->bex.addForce(id1,-f3); massSpring->bex.addForce(id2, f3); #else - static_cast<Force*> ( physicalActions->find( id1 , actionForce->getClassIndex() ).get() )->force -= f3; - static_cast<Force*> ( physicalActions->find( id2 , actionForce->getClassIndex() ).get() )->force += f3; + static_cast<Force*> ( massSpring->physicalActions->find( id1 , actionForce->getClassIndex() ).get() )->force -= f3; + static_cast<Force*> ( massSpring->physicalActions->find( id2 , actionForce->getClassIndex() ).get() )->force += f3; #endif } } _______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp _______________________________________________ yade-dev mailing list [email protected] https://lists.berlios.de/mailman/listinfo/yade-dev
