yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-04-10 10:49:36 +0200 (Fri, 10 Apr 2009)
New Revision: 1753

Modified:
   trunk/SConstruct
   trunk/extra/Brefcom.cpp
   trunk/extra/usct/UniaxialStrainControlledTest.cpp
   trunk/gui/py/yade-multi
   trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp
Log:
1. Fix race condition in brefcom (omp critical section)
2. Remove -ffast-math from optimized build, since NaNs do not work properly (e.g. isnan(NaN)==false!)



Modified: trunk/SConstruct
===================================================================
--- trunk/SConstruct	2009-04-09 10:15:52 UTC (rev 1752)
+++ trunk/SConstruct	2009-04-10 08:49:36 UTC (rev 1753)
@@ -394,7 +394,7 @@
 else: env.Append(CXXFLAGS='-O2')
 if env['openmp']: env.Append(CXXFLAGS='-fopenmp',LIBS='gomp',CPPDEFINES='YADE_OPENMP')
 if env['optimize']:
-	env.Append(CXXFLAGS=Split('-O3 -ffast-math -march=%s'%env['march']),
+	env.Append(CXXFLAGS=Split('-O3 -march=%s'%env['march']),
 		CPPDEFINES=[('YADE_CAST','static_cast'),('YADE_PTR_CAST','static_pointer_cast'),'NDEBUG'])
 	# NDEBUG is used in /usr/include/assert.h: when defined, asserts() are no-ops
 

Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp	2009-04-09 10:15:52 UTC (rev 1752)
+++ trunk/extra/Brefcom.cpp	2009-04-10 08:49:36 UTC (rev 1753)
@@ -151,21 +151,21 @@
 		Real df=(c*N*exp(N*ret)+exp(ret))/aux;
 		ret-=f/df;
 	}
-	LOG_FATAL("No convergence, c="<<c<<", ret="<<ret<<", f="<<f);
+	LOG_FATAL("No convergence after "<<maxIter<<" iters; c="<<c<<", N="<<N<<", ret="<<ret<<", f="<<f);
 	throw runtime_error("ef2_Spheres_Brefcom_BrefcomLaw::solveBeta failed to converge.");
 }
 
 Real BrefcomContact::computeDmgOverstress(Real dt){
 	if(dmgStrain>=epsN*omega){ // unloading, no viscous stress
 		dmgStrain=epsN*omega;
-		LOG_DEBUG("Elastic/unloading, no viscous overstress");
+		LOG_TRACE("Elastic/unloading, no viscous overstress");
 		return 0.;
 	}
 	Real c=epsCrackOnset*(1-omega)*pow(dmgTau/dt,dmgRateExp)*pow(epsN*omega-dmgStrain,dmgRateExp-1.);
 	Real beta=solveBeta(c,dmgRateExp);
 	Real deltaDmgStrain=(epsN*omega-dmgStrain)*exp(beta);
 	dmgStrain+=deltaDmgStrain;
-	LOG_DEBUG("deltaDmgStrain="<<deltaDmgStrain<<", viscous overstress "<<(epsN*omega-dmgStrain)*E);
+	LOG_TRACE("deltaDmgStrain="<<deltaDmgStrain<<", viscous overstress "<<(epsN*omega-dmgStrain)*E);
 	/* σN=Kn(εN-εd); dmgOverstress=σN-(1-ω)*Kn*εN=…=Kn(ω*εN-εd) */
 	return (epsN*omega-dmgStrain)*E;
 }

Modified: trunk/extra/usct/UniaxialStrainControlledTest.cpp
===================================================================
--- trunk/extra/usct/UniaxialStrainControlledTest.cpp	2009-04-09 10:15:52 UTC (rev 1752)
+++ trunk/extra/usct/UniaxialStrainControlledTest.cpp	2009-04-10 08:49:36 UTC (rev 1753)
@@ -49,7 +49,7 @@
 	assert(originalLength>0 && !isnan(originalLength));
 
 	assert(!isnan(strainRate) || !isnan(absSpeed));
-	if(strainRate==0){ strainRate=absSpeed/originalLength; LOG_INFO("COmputed new strainRate "<<strainRate); }
+	if(strainRate==0){ strainRate=absSpeed/originalLength; LOG_INFO("Computed new strainRate "<<strainRate); }
 
 	initAccelTime_s=initAccelTime>=0 ? initAccelTime : Omega::instance().getTimeStep()*(-initAccelTime);
 	LOG_INFO("Strain speed will be "<<absSpeed<<", strain rate "<<strainRate<<", will be reached after "<<initAccelTime_s<<"s ("<<initAccelTime_s/Omega::instance().getTimeStep()<<" steps).");

Modified: trunk/gui/py/yade-multi
===================================================================
--- trunk/gui/py/yade-multi	2009-04-09 10:15:52 UTC (rev 1752)
+++ trunk/gui/py/yade-multi	2009-04-10 08:49:36 UTC (rev 1753)
@@ -68,7 +68,7 @@
 parser.add_option('--nice',dest='nice',type='int',help='Nice value of spawned jobs (default: 10)',default=10)
 parser.add_option('--executable',dest='executable',help='Name of the program to run (default: %s)'%sys.argv[0][:-6],default=sys.argv[0][:-6],metavar='FILE') ## strip the '-multi' extension
 parser.add_option('--gnuplot',dest='gnuplotOut',help='Gnuplot file where gnuplot from all jobs should be put together',default=None,metavar='FILE')
-parser.add_option('--dry-run',dest='dryRun',help='Do not actually run (useful for getting gnuplot only, for instance)',default=False)
+parser.add_option('--dry-run',action='store_true',dest='dryRun',help='Do not actually run (useful for getting gnuplot only, for instance)',default=False)
 opts,args=parser.parse_args()
 logFormat,lineList,maxJobs,nice,executable,gnuplotOut,dryRun=opts.logFormat,opts.lineList,opts.maxJobs,opts.nice,opts.executable,opts.gnuplotOut,opts.dryRun
 
@@ -155,8 +155,8 @@
 				break
 	preamble,plots='',[]
 	for job in jobs:
-		if not job.plot:
-			# warn here
+		if not 'plot' in job.__dict__:
+			print "WARN: No plot found for job "+job.id
 			continue
 		for l in file(job.plot):
 			if l.startswith('plot'):

Modified: trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp	2009-04-09 10:15:52 UTC (rev 1752)
+++ trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp	2009-04-10 08:49:36 UTC (rev 1753)
@@ -179,7 +179,8 @@
 			the findX, findY, findZ take almost the totality of the time.
 			Parallelizing those is vastly beneficial (almost 3x speed increase, which can be quite sensible as the initial
 			findOverlappingBB is really slow http://yade.wikia.com/wiki/Colliders_performace and is done in 3
-			orthogonal directions. Therefore, it is enabled by default.
+			orthogonal directions. Therefore, it is enabled by default. updateOverlappingBBSet must be protected by 
+			critical section, since it is called from all threads!
 			
 			Now sortX is right before findX etc, in the same openMP section. Beware that timingDeltas will give garbage
 			results if used in parallelized code.
@@ -250,33 +251,35 @@
 
 /* Note that this function is called only for bodies that actually overlap along some axis */
 void PersistentSAPCollider::updateOverlapingBBSet(int id1,int id2){
- 	// look if the pair (id1,id2) already exists in the overlappingBB collection
-	const shared_ptr<Interaction>& interaction=transientInteractions->find(body_id_t(id1),body_id_t(id2));
-	bool found=(interaction!=0);//Bruno's Hack
-	// if there is persistent interaction, we will not create transient one!
-	
-	// test if the AABBs of the spheres number "id1" and "id2" are overlapping
-	int offset1=3*id1, offset2=3*id2;
-	const shared_ptr<Body>& b1(Body::byId(body_id_t(id1),rootBody)), b2(Body::byId(body_id_t(id2),rootBody));
-	bool overlap =
-		Collider::mayCollide(b1.get(),b2.get()) &&
-		// AABB collisions: 
-		!(
-			maxima[offset1  ]<minima[offset2  ] || maxima[offset2  ]<minima[offset1  ] || 
-			maxima[offset1+1]<minima[offset2+1] || maxima[offset2+1]<minima[offset1+1] || 
-			maxima[offset1+2]<minima[offset2+2] || maxima[offset2+2]<minima[offset1+2]);
-	// inserts the pair p=(id1,id2) if the two AABB overlaps and if p does not exists in the overlappingBB
-	//if((id1==0 && id2==1) || (id1==1 && id2==0)) LOG_DEBUG("Processing #0 #1");
-	//if(interaction&&!interaction->isReal){ LOG_DEBUG("Unreal interaction #"<<id1<<"=#"<<id2<<" (overlap="<<overlap<<", haveDistantTransient="<<haveDistantTransient<<")");}
-	if(overlap && !found){
-		//LOG_DEBUG("Creating interaction #"<<id1<<"=#"<<id2);
-		transientInteractions->insert(body_id_t(id1),body_id_t(id2));
+	#pragma omp critical
+	{
+		// look if the pair (id1,id2) already exists in the overlappingBB collection
+		const shared_ptr<Interaction>& interaction=transientInteractions->find(body_id_t(id1),body_id_t(id2));
+		bool found=(bool)interaction;
+		
+		// test if the AABBs of the spheres number "id1" and "id2" are overlapping
+		int offset1=3*id1, offset2=3*id2;
+		const shared_ptr<Body>& b1(Body::byId(body_id_t(id1),rootBody)), b2(Body::byId(body_id_t(id2),rootBody));
+		bool overlap =
+			Collider::mayCollide(b1.get(),b2.get()) &&
+			// AABB collisions: 
+			!(
+				maxima[offset1  ]<minima[offset2  ] || maxima[offset2  ]<minima[offset1  ] || 
+				maxima[offset1+1]<minima[offset2+1] || maxima[offset2+1]<minima[offset1+1] || 
+				maxima[offset1+2]<minima[offset2+2] || maxima[offset2+2]<minima[offset1+2]);
+		// inserts the pair p=(id1,id2) if the two AABB overlaps and if p does not exists in the overlappingBB
+		//if((id1==0 && id2==1) || (id1==1 && id2==0)) LOG_DEBUG("Processing #0 #1");
+		//if(interaction&&!interaction->isReal){ LOG_DEBUG("Unreal interaction #"<<id1<<"=#"<<id2<<" (overlap="<<overlap<<", haveDistantTransient="<<haveDistantTransient<<")");}
+		if(overlap && !found){
+			//LOG_DEBUG("Creating interaction #"<<id1<<"=#"<<id2);
+			transientInteractions->insert(body_id_t(id1),body_id_t(id2));
+		}
+		// removes the pair p=(id1,id2) if the two AABB do not overlapp any more and if p already exists in the overlappingBB
+		else if(!overlap && found && (haveDistantTransient ? !interaction->isReal : true) ){
+			//LOG_DEBUG("Erasing interaction #"<<id1<<"=#"<<id2<<" (isReal="<<interaction->isReal<<")");
+			transientInteractions->erase(body_id_t(id1),body_id_t(id2));
+		}
 	}
-	// removes the pair p=(id1,id2) if the two AABB do not overlapp any more and if p already exists in the overlappingBB
-	else if(!overlap && found && (haveDistantTransient ? !interaction->isReal : true) ){
-		//LOG_DEBUG("Erasing interaction #"<<id1<<"=#"<<id2<<" (isReal="<<interaction->isReal<<")");
-		transientInteractions->erase(body_id_t(id1),body_id_t(id2));
-	}
 }