[cig-commits] r14244 - long/3D/SNAC/trunk/Snac/plugins/hillSlope

cstark at geodynamics.org cstark at geodynamics.org
Sat Mar 7 13:40:53 PST 2009


Author: cstark
Date: 2009-03-07 13:40:52 -0800 (Sat, 07 Mar 2009)
New Revision: 14244

Modified:
   long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c
   long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h
   long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.c
   long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.c
Log:
Tracking, dumping improved, extended

1) Cleaned handling of elastic stability (eqm) tracking
2) Added post-eqm dump frequency variable
3) Ensured dump of elastic eqm state if solving that only

Details:

1a) Tidied flagging and flag handling in Track.c
1b) Ensured plastic weak points assignments are not carried out in CreateWeakPoints.c if solving elastic eqm only (bug fix)

2) Added new contextExt variable plasticDeformationDumpFreq that sets rate of state dumps once elastic eqm has been reached.

3a) Track.c now changes dump frequency when elastic eqm reached.  An MPI call ensures all threads have a consensus value of contextExt->dumpEvery once this change is effected.
3b) TO DO: need to alter this behavior a bit to have dump freq set to one if solving elastic eqm only.




Modified: long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c	2009-03-06 21:12:23 UTC (rev 14243)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c	2009-03-07 21:40:52 UTC (rev 14244)
@@ -83,6 +83,9 @@
 	contextExt->rngSeed = Dictionary_Entry_Value_AsUnsignedInt(
 		Dictionary_GetDefault( context->dictionary, "rngSeed", 
 				       Dictionary_Entry_Value_FromUnsignedInt( 1 ) ) );
+	contextExt->plasticDeformationDumpFreq = Dictionary_Entry_Value_AsUnsignedInt(
+		Dictionary_GetDefault( context->dictionary, "plasticDeformationDumpFreq", 
+				       Dictionary_Entry_Value_FromUnsignedInt( 100 ) ) );
 
 	contextExt->fractionWeakPoints = Dictionary_Entry_Value_AsDouble(
 		Dictionary_GetDefault( context->dictionary, "fractionWeakPoints", 

Modified: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h	2009-03-06 21:12:23 UTC (rev 14243)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h	2009-03-07 21:40:52 UTC (rev 14244)
@@ -48,6 +48,7 @@
 		/* For _SnacHillSlope_Top2BottomSweep condition function */
 		double				slopeAngle;
 		unsigned int			rngSeed;
+		unsigned int			plasticDeformationDumpFreq;
 		double				fractionWeakPoints;
 		double				resolveDepth;
 		double				leftFlatFraction;

Modified: long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.c	2009-03-06 21:12:23 UTC (rev 14243)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.c	2009-03-07 21:40:52 UTC (rev 14244)
@@ -127,8 +127,11 @@
 
     /*
      *  Bail now if initial elastic equilibrium has not been reached on all threads
+     *
+     *  Also bail if we're solving elastic eqm only
      */
-    if(!contextExt->consensusElasticStabilizedFlag || contextExt->seedingCompletedFlag) return;
+    if(!contextExt->consensusElasticStabilizedFlag || contextExt->seedingCompletedFlag
+       || (contextExt->solveElasticEqmOnlyFlag)) return;
 
     //    fprintf(stderr, "CWP\n");
 
@@ -289,8 +292,8 @@
 	    /*  Force each of 5*2 tetrahedra to have extra plastic strain to impose lower cohesion indirectly */
 	    for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
 		plasticElement->plasticStrain[tetra_I] =  PlasticStrainFromCohesion(material,(double)contextExt->triggerPointCohesion);
+#ifdef DEBUG
 		if(tetra_I==0) Journal_Printf(context->snacInfo, "timeStep=%d ijkt=%d %d %d %d  setting plasticE=%e\n", context->timeStep,ijk[0],ijk[1],ijk[2], tetra_I,  plasticElement->plasticStrain[tetra_I] );
-#ifdef DEBUG
 #endif
 	    }
 	}
@@ -371,10 +374,11 @@
  */
 
 void 
-ShellSort(vecPtr,idxPtr,n)
-	unsigned int vecPtr[];              /* Input unsorted data vector (e.g. random numbers) */  /* was double */
-	unsigned int idxPtr[];              /* Input index vector, to be sorted according to vecPtr */
-	unsigned long n;          	    /* Data length */
+ShellSort(
+	unsigned int vecPtr[],              /* Input unsorted data vector (e.g. random numbers) */  /* was double */
+	unsigned int idxPtr[],              /* Input index vector, to be sorted according to vecPtr */
+	unsigned long n          	    /* Data length */
+	)
 {
     unsigned long incr;
     incr=1;

Modified: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.c	2009-03-06 21:12:23 UTC (rev 14243)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.c	2009-03-07 21:40:52 UTC (rev 14244)
@@ -44,220 +44,225 @@
 #define FALSE 0
 #endif
 
+#define DEBUG
 #define DEBUG2
 
 void SnacHillSlope_Track( void* _context ) {
-	Snac_Context			*context = (Snac_Context*)_context;
-	SnacHillSlope_Context		*contextExt = ExtensionManager_Get(
-							context->extensionMgr,
-							context,
-							SnacHillSlope_ContextHandle );
-	Mesh				*mesh = context->mesh;
-	MeshLayout			*layout = (MeshLayout*)mesh->layout;
-	HexaMD				*decomp = (HexaMD*)layout->decomp;
-	BlockGeometry			*geometry = (BlockGeometry*)layout->elementLayout->geometry;
-	const int			full_I_node_range=decomp->nodeGlobal3DCounts[0];
-	const int			full_J_node_range=decomp->nodeGlobal3DCounts[1];
-	const int			full_K_node_range=decomp->nodeGlobal3DCounts[2];
-	int                   	        index_I,index_J,index_K;
+    Snac_Context		*context = (Snac_Context*)_context;
+    SnacHillSlope_Context	*contextExt = ExtensionManager_Get(context->extensionMgr,
+								   context,
+								   SnacHillSlope_ContextHandle );
+    Mesh			*mesh = context->mesh;
+    MeshLayout			*layout = (MeshLayout*)mesh->layout;
+    HexaMD			*decomp = (HexaMD*)layout->decomp;
+    BlockGeometry		*geometry = (BlockGeometry*)layout->elementLayout->geometry;
+    const int			full_I_node_range=decomp->nodeGlobal3DCounts[0];
+    const int			full_J_node_range=decomp->nodeGlobal3DCounts[1];
+    const int			full_K_node_range=decomp->nodeGlobal3DCounts[2];
+    int				index_I,index_J,index_K;
 
-	static double			*yGridOldPtr, *yGridOlderPtr;
-	double				max_yVelocity, max_yAcceln;
-	static double			unit_yVelocity=0.0, unit_yAcceln=0.0;
-	static double			old_max_yVelocity=0.0, old_max_yAcceln=0.0;
-/* 	static double			min_yVelocity=0.0, min_yAcceln=0.0; */
-	static char			fallingFlag=FALSE;
-	unsigned char			reachesTopFlag=FALSE;
+    static double		*yGridOldPtr, *yGridOlderPtr;
+    double			max_yVelocity, max_yAcceln;
+    static double		unit_yVelocity=0.0, unit_yAcceln=0.0;
+    static double		old_max_yVelocity=0.0, old_max_yAcceln=0.0;
+    /* 	static double			min_yVelocity=0.0, min_yAcceln=0.0; */
+    static char			fallingFlag=FALSE, doneTrackingFlag=FALSE;
+    char			reachesTopFlag=FALSE;
 
-	const double			trackLevel=(double)contextExt->trackLevel;
-	const double			startThreshold=(contextExt->startThreshold>=0.0 ? contextExt->startThreshold : 1e-2);
-	const double			stopThreshold=(contextExt->stopThreshold>=0.0 ? contextExt->stopThreshold : 1e-3);
-	char eflag, ceflag;
-	
+    const double		trackLevel=(double)contextExt->trackLevel;
+    const double		startThreshold=(contextExt->startThreshold>=0.0 ? contextExt->startThreshold : 1e-2);
+    const double		stopThreshold=(contextExt->stopThreshold>=0.0 ? contextExt->stopThreshold : 1e-3);
+    int				maxTimeSteps=context->maxTimeSteps, dumpEvery=context->dumpEvery;
 
-/* 	if (context->timeStep % context->dumpEvery == 0) { */
-/* 	    Journal_Printf( context->snacInfo,"timeStep=%d (in track)\n", context->timeStep ); */
-/* 	} */
+    /* 	if (context->timeStep % context->dumpEvery == 0) { */
+    /* 	    Journal_Printf( context->snacInfo,"timeStep=%d (in track)\n", context->timeStep ); */
+    /* 	} */
 
-	/*
-	 *  Bail now if all threads have reached elastic equilibrium
-	 */
-	if(contextExt->consensusElasticStabilizedFlag){
+    /*
+     *  Bail now if all threads have reached elastic equilibrium
+     */
+    if(contextExt->consensusElasticStabilizedFlag || doneTrackingFlag){
 #ifdef DEBUG
-	    fprintf(stderr,"c=%d,t=%d/%d: Consensus eqm... bailing at top of Track.c\n",context->rank, 
-		    context->timeStep,context->maxTimeSteps);
+	fprintf(stderr,"c=%d,t=%d/%d: Consensus eqm... bailing at top of Track.c\n",context->rank, 
+		context->timeStep,context->maxTimeSteps);
 #endif
-	    return;
-	}
+	return;
+    }
 
-	/*
-	 *  Set up a tracking grids (slices of mesh) to allow t instance to be compared with t-1, t-2 instances
-	 */
-	if(context->timeStep==1) {
-	    yGridOldPtr=(double *)malloc((size_t)(full_I_node_range*full_K_node_range*sizeof(double)));
-	    yGridOlderPtr=(double *)malloc((size_t)(full_I_node_range*full_K_node_range*sizeof(double)));
-	    for(index_I = 0; index_I < full_I_node_range; index_I++) {
-		for(index_K = 0; index_K < full_K_node_range; index_K++) {
-		    *(yGridOlderPtr+index_I + full_I_node_range*index_K)=-1.0E19;
-		    *(yGridOldPtr+index_I + full_I_node_range*index_K)=-1.0E19;
-		}
+    /*
+     *  Set up a tracking grids (slices of mesh) to allow t instance to be compared with t-1, t-2 instances
+     */
+    if(context->timeStep==1) {
+	yGridOldPtr=(double *)malloc((size_t)(full_I_node_range*full_K_node_range*sizeof(double)));
+	yGridOlderPtr=(double *)malloc((size_t)(full_I_node_range*full_K_node_range*sizeof(double)));
+	for(index_I = 0; index_I < full_I_node_range; index_I++) {
+	    for(index_K = 0; index_K < full_K_node_range; index_K++) {
+		*(yGridOlderPtr+index_I + full_I_node_range*index_K)=-1.0E19;
+		*(yGridOldPtr+index_I + full_I_node_range*index_K)=-1.0E19;
 	    }
 	}
+    }
 
-	/*
-	 * Compare t instance with t-1, t-2 instances, calculating slice vel, acceln and identifying maxima
-	 */
-	max_yVelocity=0.0;
-	max_yAcceln=0.0;
-	for(index_I = 0; index_I < full_I_node_range; index_I++) {
-	    for(index_K = 0; index_K < full_K_node_range; index_K++) {
-		Node_GlobalIndex		node_gI;
-		Node_LocalIndex			node_lI;
-		Coord				*coordPtr = 0;
-		double				node_yElevation, node_yVelocity, node_yAcceln;
-		double				*tmp_yGridOldPtr, *tmp_yGridOlderPtr;
+    /*
+     * Compare t instance with t-1, t-2 instances, calculating slice vel, acceln and identifying maxima
+     */
+    max_yVelocity=0.0;
+    max_yAcceln=0.0;
+    for(index_I = 0; index_I < full_I_node_range; index_I++) {
+	for(index_K = 0; index_K < full_K_node_range; index_K++) {
+	    Node_GlobalIndex		node_gI;
+	    Node_LocalIndex			node_lI;
+	    Coord				*coordPtr = 0;
+	    double				node_yElevation, node_yVelocity, node_yAcceln;
+	    double				*tmp_yGridOldPtr, *tmp_yGridOlderPtr;
 
-		index_J=(int)(((double)full_J_node_range-1.0)*(1.0-trackLevel));
-		node_gI = index_I + full_I_node_range*index_J + full_I_node_range*full_J_node_range*index_K;
-		node_lI = Mesh_NodeMapGlobalToLocal( mesh, node_gI );
+	    index_J=(int)(((double)full_J_node_range-1.0)*(1.0-trackLevel));
+	    node_gI = index_I + full_I_node_range*index_J + full_I_node_range*full_J_node_range*index_K;
+	    node_lI = Mesh_NodeMapGlobalToLocal( mesh, node_gI );
 
+	    /* 
+	     * If a local node, read its elevation, if not, give a dummy value 
+	     */
+	    if( node_lI < context->mesh->nodeLocalCount ) { 
 		/* 
-		 * If a local node, read its elevation, if not, give a dummy value 
+		 *  A local node 
 		 */
-		if( node_lI < context->mesh->nodeLocalCount ) { 
-		    /* 
-		     *  A local node 
-		     */
-		    coordPtr = Snac_NodeCoord_P( context, node_lI );
-		    node_yElevation = (*coordPtr)[1]+(*coordPtr)[1]+(*coordPtr)[1];
-		    reachesTopFlag |= TRUE;
-		}
-		else {
-		    /*
-		     *  Not a local node - handled by another CPU thread
-		     */
-		    node_yElevation=-1.0E29;
-		}
-		tmp_yGridOldPtr = yGridOldPtr+index_I + full_I_node_range*index_K;
-		tmp_yGridOlderPtr = yGridOlderPtr+index_I + full_I_node_range*index_K;
+		coordPtr = Snac_NodeCoord_P( context, node_lI );
+		node_yElevation = (*coordPtr)[1]+(*coordPtr)[1]+(*coordPtr)[1];
+		reachesTopFlag |= TRUE;
+	    }
+	    else {
 		/*
-		 * Start assessing the vertical surface motions after two time steps
+		 *  Not a local node - handled by another CPU thread
 		 */
-		node_yVelocity = node_yElevation-*tmp_yGridOldPtr;
-		node_yAcceln = node_yVelocity-(*tmp_yGridOldPtr-*tmp_yGridOlderPtr);
-		if(context->timeStep>=3){
-		    if(fabs(node_yVelocity)>max_yVelocity)
-			max_yVelocity = fabs(node_yVelocity);
-		    if(fabs(node_yAcceln)>max_yAcceln)
-			max_yAcceln = fabs(node_yAcceln);
-		}
-		/*
-		 * Record this elevation field and push previous back to "OlderPtr" array
-		 */
-		*tmp_yGridOlderPtr = *tmp_yGridOldPtr;
-		*tmp_yGridOldPtr = node_yElevation;
+		node_yElevation=-1.0E29;
 	    }
+	    tmp_yGridOldPtr = yGridOldPtr+index_I + full_I_node_range*index_K;
+	    tmp_yGridOlderPtr = yGridOlderPtr+index_I + full_I_node_range*index_K;
+	    /*
+	     * Start assessing the vertical surface motions after two time steps
+	     */
+	    node_yVelocity = node_yElevation-*tmp_yGridOldPtr;
+	    node_yAcceln = node_yVelocity-(*tmp_yGridOldPtr-*tmp_yGridOlderPtr);
+	    if(context->timeStep>=3){
+		if(fabs(node_yVelocity)>max_yVelocity)
+		    max_yVelocity = fabs(node_yVelocity);
+		if(fabs(node_yAcceln)>max_yAcceln)
+		    max_yAcceln = fabs(node_yAcceln);
+	    }
+	    /*
+	     * Record this elevation field and push previous back to "OlderPtr" array
+	     */
+	    *tmp_yGridOlderPtr = *tmp_yGridOldPtr;
+	    *tmp_yGridOldPtr = node_yElevation;
 	}
+    }
 
 #ifdef DEBUG2
-	    fprintf(stderr,
-		    "c=%d,t=%d/%d:  reachesTop=%d consensusElasticStabilized=%d  elasticStabilized=%d   startedTracking=%d:  max_vel=%g  unit_vel=%g\n",
-		    context->rank, context->timeStep, context->maxTimeSteps, reachesTopFlag, 
-		    contextExt->consensusElasticStabilizedFlag, contextExt->elasticStabilizedFlag, contextExt->startedTrackingFlag, 
-		    max_yVelocity, unit_yVelocity ); 
+    fprintf(stderr,
+	    "c=%d,t=%d/%d:  reachesTop=%d consensusElasticStabilized=%d  elasticStabilized=%d   startedTracking=%d:  max_vel=%g  unit_vel=%g\n",
+	    context->rank, context->timeStep, context->maxTimeSteps, reachesTopFlag, 
+	    contextExt->consensusElasticStabilizedFlag, contextExt->elasticStabilizedFlag, contextExt->startedTrackingFlag, 
+	    max_yVelocity, unit_yVelocity ); 
 #endif
 
+    /*
+     *  If none of the top-surface nodes are local, set flags to pretend we're equilibrated on this thread
+     */
+    if(!reachesTopFlag) {
+	contextExt->startedTrackingFlag = TRUE;
+	contextExt->elasticStabilizedFlag = TRUE;
+#ifdef DEBUG
+	fprintf(stderr,"c=%d,t=%d/%d:  Doesn't reach top\n",context->rank, context->timeStep, context->maxTimeSteps);
+#endif
+    } else {
 	/*
-	 *  If none of the top-surface nodes are local, set flags to pretend we're equilibrated on this thread
+	 * Now deprecated: estimate unit rates of motion for later comparison with falling rates
 	 */
-	if(!reachesTopFlag) {
-	    contextExt->startedTrackingFlag = TRUE;
-	    contextExt->elasticStabilizedFlag = TRUE;
-#ifdef DEBUG
-	    fprintf(stderr,"c=%d,t=%d/%d:  Doesn't reach top\n",context->rank, context->timeStep, context->maxTimeSteps);
-#endif
-	} else {
-	    /*
-	     * Now deprecated: estimate unit rates of motion for later comparison with falling rates
-	     */
-	    if(unit_yVelocity==0.0 && max_yVelocity>0.0)
-		unit_yVelocity = max_yVelocity;		
-	    if(unit_yAcceln==0.0 && max_yAcceln>0.0)
-		unit_yAcceln = max_yAcceln;	
-	    if(!contextExt->startedTrackingFlag && max_yVelocity>=unit_yVelocity && context->timeStep>=4) 
-		contextExt->startedTrackingFlag=TRUE;
+	if(unit_yVelocity==0.0 && max_yVelocity>0.0)
+	    unit_yVelocity = max_yVelocity;		
+	if(unit_yAcceln==0.0 && max_yAcceln>0.0)
+	    unit_yAcceln = max_yAcceln;	
+	if(!contextExt->startedTrackingFlag && max_yVelocity>=unit_yVelocity && context->timeStep>=4) 
+	    contextExt->startedTrackingFlag=TRUE;
 
 #ifdef DEBUG
-	    fprintf(stderr,"c=%d,t=%d/%d: Does reach top - check if equilibrating\n",context->rank, context->timeStep, context->maxTimeSteps);
+	fprintf(stderr,"c=%d,t=%d/%d: Does reach top - check if equilibrating\n",context->rank, context->timeStep, context->maxTimeSteps);
 #endif
-	    /*
-	     *  If surface change is slowing and slowly enough, flag that elastic eqm has been reached on this thread
-	     */
-	    if(contextExt->startedTrackingFlag){
-		fallingFlag = CheckFallingFn(max_yVelocity,max_yAcceln,old_max_yVelocity,old_max_yAcceln);
-		/* 	    if(CheckStabilizingFn(max_yVelocity/unit_yVelocity, max_yAcceln/unit_yAcceln, stopThreshold, fallingFlag)==TRUE */
-		if( !contextExt->elasticStabilizedFlag
-		    && CheckStabilizingFn(max_yVelocity, max_yAcceln, stopThreshold, fallingFlag)==TRUE
-		    && context->maxTimeSteps!=context->timeStep ) {
-		    /*
-		     *  Stabilizing on this thread
-		     */
+	/*
+	 *  If surface change is slowing and slowly enough, flag that elastic eqm has been reached on this thread
+	 */
+	if(contextExt->startedTrackingFlag){
+	    fallingFlag = CheckFallingFn(max_yVelocity,max_yAcceln,old_max_yVelocity,old_max_yAcceln);
+	    /* 	    if(CheckStabilizingFn(max_yVelocity/unit_yVelocity, max_yAcceln/unit_yAcceln, stopThreshold, fallingFlag)==TRUE */
+	    if( !contextExt->elasticStabilizedFlag
+		&& CheckStabilizingFn(max_yVelocity, max_yAcceln, stopThreshold, fallingFlag)==TRUE
+		&& context->maxTimeSteps!=context->timeStep ) {
+		/*
+		 *  Stabilizing on this thread
+		 */
 #ifdef DEBUG
-		  fprintf(stderr,"c=%d,t=%d/%d: Locally report we're equilibrating\n",context->rank, context->timeStep, context->maxTimeSteps);
+		fprintf(stderr,"c=%d,t=%d/%d: Locally report we're equilibrating\n",context->rank, context->timeStep, context->maxTimeSteps);
 #endif
-		    contextExt->elasticStabilizedFlag = TRUE;
-		}
+		contextExt->elasticStabilizedFlag = TRUE;
 	    }
 	}
-	/*
-	 *  Decide whether to stop or to continue simulation
-	 */
+    }
+    /*
+     *  Decide whether to stop or to continue simulation
+     */
 #ifdef DEBUG
-	fprintf(stderr,"c=%d,t=%d/%d: Checking consensus... %d\n",context->rank, context->timeStep, context->maxTimeSteps,
-		contextExt->consensusElasticStabilizedFlag);
+    fprintf(stderr,"c=%d,t=%d/%d: Checking consensus... %d\n",context->rank, context->timeStep, context->maxTimeSteps,
+	    contextExt->consensusElasticStabilizedFlag);
 #endif
+    /*
+     *  Check all threads to see if global equilibration has been reached
+     */
+    MPI_Allreduce( &(contextExt->elasticStabilizedFlag), &(contextExt->consensusElasticStabilizedFlag), 
+		   1, MPI_CHAR, MPI_MIN, context->communicator );
+#ifdef DEBUG
+    fprintf(stderr,"c=%d,t=%d/%d: ... revised consensus=%d\n",context->rank, context->timeStep, context->maxTimeSteps,
+	    contextExt->consensusElasticStabilizedFlag);
+#endif
+    /* 	if(contextExt->startedTrackingFlag  && !doneTrackingFlag ){ */
+    /*
+     *  If all threads agree to elastic eqm, and we only want to run to this point, 
+     *    tell the simulation to stop at this time step
+     */
+    if(contextExt->consensusElasticStabilizedFlag) {
 	/*
-	 *  Check all threads to see if global equilibration has been reached
+	 * Stop the simulation by bringing forward the max time steps to this step.
+	 * In addition, force a dump of this model state by changing dump freq to 1.
 	 */
-	//	    MPI_Allreduce( &contextExt->elasticStabilizedFlag, &contextExt->consensusElasticStabilizedFlag, 
-	//	   1, MPI_INT, MPI_LAND, context->communicator );
-	eflag=contextExt->elasticStabilizedFlag;
-	MPI_Allreduce( &eflag, &ceflag, 
-		       1, MPI_CHAR, MPI_MIN, context->communicator );
+	dumpEvery=contextExt->plasticDeformationDumpFreq;
+	if(contextExt->solveElasticEqmOnlyFlag) {
+	    maxTimeSteps=context->timeStep+1;
+	    doneTrackingFlag=TRUE;
+	}
 #ifdef DEBUG
-	fprintf(stderr,"c=%d,t=%d/%d: ... revised consensus= %d->%d\n",context->rank, context->timeStep, context->maxTimeSteps,
-		ceflag,
-		contextExt->consensusElasticStabilizedFlag);
+	fprintf(stderr,"c=%d, t=%d/%d: Stopping run at t= %d\n",context->rank, context->timeStep, 
+		maxTimeSteps,context->maxTimeSteps);
 #endif
-	if(contextExt->startedTrackingFlag){
-	    /*
-	     *  If all threads agree to elastic eqm, and we only want to run to this point, 
-	     *    tell the simulation to stop at this time step
-	     */
-	  contextExt->consensusElasticStabilizedFlag=ceflag;
-#ifdef DEBUG
-	fprintf(stderr,"c=%d,t=%d/%d: ... revised consensus= %d->%d\n",context->rank, context->timeStep, context->maxTimeSteps,
-		ceflag,
-		contextExt->consensusElasticStabilizedFlag);
-#endif
+    }
+    /* 	} */
+    /*
+     *  Ensure that all threads are in agreement
+     */
+    if(contextExt->consensusElasticStabilizedFlag) {
+	MPI_Allreduce( &(maxTimeSteps), &(context->maxTimeSteps), 1, MPI_INT, MPI_MIN, context->communicator );
+	MPI_Allreduce( &(dumpEvery), &(context->dumpEvery), 1, MPI_INT, MPI_MIN, context->communicator );
+    }
 
-	  if(contextExt->consensusElasticStabilizedFlag) {
-	    if(contextExt->solveElasticEqmOnlyFlag) {
 #ifdef DEBUG
-	      fprintf(stderr,"c=%d, t=%d/%d: Stopping run at t= %d\n",context->rank, context->timeStep, 
-		      context->maxTimeSteps,context->maxTimeSteps);
+    fprintf(stderr,"c=%d, t=%d/%d: Exiting Track.c with dumpFreq=%d\n",
+	    context->rank, context->timeStep, context->maxTimeSteps, context->dumpEvery);
 #endif
-	      context->maxTimeSteps=context->timeStep/*+context->dumpEvery*/;
-	    }
-	  }
-	}
+    /*
+     * Record the current mesh slice velocity and acceln for use next iteration
+     */
+    old_max_yVelocity = max_yVelocity;
+    old_max_yAcceln = max_yAcceln;
 
-	/*
-	 * Record the current mesh slice velocity and acceln for use next iteration
-	 */
-	old_max_yVelocity = max_yVelocity;
-	old_max_yAcceln = max_yAcceln;
-
 }
 
 



More information about the CIG-COMMITS mailing list