[cig-commits] r14656 - in long/3D/SNAC/trunk/Snac: libSnac/src plugins/hillSlope plugins/maxwell plugins/plastic plugins/remesher plugins/restarter plugins/restarter_old plugins/temperature plugins/viscoplastic plugins/winkler snac2vtk
echoi at geodynamics.org
echoi at geodynamics.org
Fri Apr 10 13:43:46 PDT 2009
Author: echoi
Date: 2009-04-10 13:43:44 -0700 (Fri, 10 Apr 2009)
New Revision: 14656
Modified:
long/3D/SNAC/trunk/Snac/libSnac/src/Context.c
long/3D/SNAC/trunk/Snac/libSnac/src/Context.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.c
long/3D/SNAC/trunk/Snac/plugins/maxwell/ConstructExtensions.c
long/3D/SNAC/trunk/Snac/plugins/maxwell/Context.h
long/3D/SNAC/trunk/Snac/plugins/maxwell/Output.c
long/3D/SNAC/trunk/Snac/plugins/maxwell/Output.h
long/3D/SNAC/trunk/Snac/plugins/maxwell/Register.c
long/3D/SNAC/trunk/Snac/plugins/plastic/Constitutive.c
long/3D/SNAC/trunk/Snac/plugins/plastic/ConstructExtensions.c
long/3D/SNAC/trunk/Snac/plugins/plastic/Context.h
long/3D/SNAC/trunk/Snac/plugins/plastic/Output.c
long/3D/SNAC/trunk/Snac/plugins/plastic/Output.h
long/3D/SNAC/trunk/Snac/plugins/plastic/Register.c
long/3D/SNAC/trunk/Snac/plugins/plastic/Remesh.c
long/3D/SNAC/trunk/Snac/plugins/plastic/makefile
long/3D/SNAC/trunk/Snac/plugins/remesher/InitialConditions.c
long/3D/SNAC/trunk/Snac/plugins/restarter/InitialConditions.c
long/3D/SNAC/trunk/Snac/plugins/restarter_old/InitialConditions.c
long/3D/SNAC/trunk/Snac/plugins/temperature/ConstructExtensions.c
long/3D/SNAC/trunk/Snac/plugins/temperature/Context.h
long/3D/SNAC/trunk/Snac/plugins/temperature/Output.c
long/3D/SNAC/trunk/Snac/plugins/temperature/Output.h
long/3D/SNAC/trunk/Snac/plugins/temperature/Register.c
long/3D/SNAC/trunk/Snac/plugins/temperature/VariableConditions.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Constitutive.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic/ConstructExtensions.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Context.h
long/3D/SNAC/trunk/Snac/plugins/viscoplastic/InitialConditions.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Output.c
long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Output.h
long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Register.c
long/3D/SNAC/trunk/Snac/plugins/winkler/Force.c
long/3D/SNAC/trunk/Snac/plugins/winkler/InitialConditions.c
long/3D/SNAC/trunk/Snac/plugins/winkler/Output.c
long/3D/SNAC/trunk/Snac/snac2vtk/snac2vtk.c
Log:
Dumping and checkpointing have been separated.
* Dumping=writing outputs in the usual sense.
Checkpointing=Saving data for restarting.
* Related input parameters:
- dumpEvery: the frequency of dumping.
- checkpointEvery: the frequency of checkpointing.
-- 0 by default.
-- Even when it's 0, the last step's data are *always* checkpointed.
* The repeated theme of all the committed changes is
- *_Dump* functions are replace with *_Write*.
- Within each *_Write* function,
*_Dump* is called if isTimeToDump() is True;
*_Checkpoint* is called if isTimeToCheckpoint() is True.
- This way, the data saving frequencies can be centrally
controlled by these two functions
(i.e., isTimeTo* defined in libSnac/src/Context.c).
* When dumping, stress tensors are averaged over each element before written to files. The file size has been greatly reduced as a result.
When checkpointing, stress vectors are still recorded for each tetrahedra.
- snac2vtk has been modified to be compatible with this change.
* Restarting has yet to be worked on.
- As a first step, "restartStep" has been replaced with "restartTimestep". The latter is defined on the AbstractContext level.
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Context.c
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Context.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Context.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -297,7 +297,7 @@
/* Snac_Context info */
/* Although we set this for safely (use of it will give 0s), it shouldn't be used... the new value should be set before
hand */
- self->restartStep = Dictionary_Entry_Value_AsUnsignedInt(
+ self->restartTimestep = Dictionary_Entry_Value_AsUnsignedInt(
Dictionary_GetDefault( self->dictionary, "restartStep", Dictionary_Entry_Value_FromUnsignedInt( 0 ) ) );
tmpStr = Dictionary_Entry_Value_AsString( Dictionary_Get( self->dictionary, "dtType" ) );
@@ -400,6 +400,15 @@
self->snacError,
"\"%s\" failed to open file for writing", tmpBuf );
}
+ if( self->checkpointEvery > 0 ) {
+ sprintf( tmpBuf, "%s/checkpointTimeStep.%u", self->outputPath, self->rank );
+ if( (self->checkpointTimeStepInfo = fopen( tmpBuf, "w+" )) == NULL ) {
+ Journal_Firewall(
+ (ArithPointer)self->checkpointTimeStepInfo,
+ self->snacError,
+ "\"%s\" failed to open file for writing", tmpBuf );
+ }
+ }
/* Add hooks to exisiting entry points */
EntryPoint_Append(
@@ -915,17 +924,12 @@
Snac_Context* self = (Snac_Context*)context;
Element_LocalIndex element_lI;
- /* if restarting, reset the time step counter */
- if( self->restartStep > 0 )
- self->timeStep = self->restartStep;
-
#ifdef DEBUG
- fprintf(stderr, "TimeStepZero: restartStep=%d, timeStep=%d\n", self->restartStep, self->timeStep);
+ fprintf(stderr, "TimeStepZero: restartTimestep=%d, timeStep=%d\n", self->restartTimestep, self->timeStep);
#endif
- _Snac_Context_InitDump( self );
+ _Snac_Context_InitOutput( self );
-
if( self->rank == 0 ) Journal_Printf( self->snacInfo, "In: %s\n", __func__ );
if( self->rank == 0 ) Journal_Printf(
self->info,
@@ -934,10 +938,10 @@
self->currentTime );
/* Write out timeStep info */
- DumpLoopInfo( self );
+ _Snac_Context_WriteLoopInfo( self );
- if( self->restartStep == 0 ) {
- for( element_lI = 0; element_lI < self->mesh->elementLocalCount && self->restartStep == 0; element_lI++ ) {
+ if( self->restartTimestep == 0 ) {
+ for( element_lI = 0; element_lI < self->mesh->elementLocalCount && self->restartTimestep == 0; element_lI++ ) {
Snac_Element* element = Snac_Element_At( self, element_lI );
Material_Index material_I = element->material_I;
Snac_Material* material = &self->materialProperty[material_I];
@@ -961,26 +965,11 @@
}
}
- Journal_Dump( self->strainRateOut, &(self->timeStep) );
-
- Journal_Dump( self->stressOut, &(self->timeStep) );
-
- Journal_Dump( self->hydroPressureOut, &(self->timeStep) );
-
- _Snac_Context_DumpStressTensor( self );
-
-
/* Update all the elements, and in the process work out this processor's minLengthScale */
KeyCall( self, self->loopElementsMomentumK, EntryPoint_VoidPtr_CallCast* )( KeyHandle(self,self->loopElementsMomentumK), self );
- Journal_Dump( self->coordOut, &(self->timeStep) );
+ _Snac_Context_WriteOutput( self );
- Journal_Dump( self->velOut, &(self->timeStep) );
-
- Journal_Dump( self->forceOut, &(self->timeStep) );
-
- Journal_Dump( self->phaseIndexOut, &(self->timeStep) );
-
KeyCall( self, self->syncK, EntryPoint_Class_VoidPtr_CallCast* )( KeyHandle(self,self->syncK), self );
/* _Snac_Context_Sync( self ); */
}
@@ -991,10 +980,6 @@
if( self->rank == 0 ) Journal_Printf( self->debug, "In: %s\n", __func__ );
- /* if restarting, reset the time step counter */
- if( self->restartStep > 0 && self->timeStep <= 1)
- self->timeStep = self->restartStep+1;
-
return self->dt;
}
@@ -1019,9 +1004,7 @@
if( self->rank == 0 ) Journal_DPrintf( self->debug, "In: %s\n", __func__ );
- if( self->timeStep == 0 || ((self->timeStep-1) % self->dumpEvery == 0) ) {
- DumpLoopInfo( self );
- }
+ _Snac_Context_WriteLoopInfo( self );
/* Perform the Snac solve loop */
/* update the energy solve */
@@ -1038,27 +1021,20 @@
Mesh_Sync( self->mesh );
}
- Journal_Dump( self->strainRateOut, &(self->timeStep) );
- Journal_Dump( self->stressOut, &(self->timeStep) );
- Journal_Dump( self->hydroPressureOut, &(self->timeStep) );
- _Snac_Context_DumpStressTensor( self );
-
KeyCall( self, self->loopNodesMomentumK, EntryPoint_VoidPtr_CallCast* )( KeyHandle(self,self->loopNodesMomentumK), self );
- Journal_Dump( self->coordOut, &(self->timeStep) );
- Journal_Dump( self->velOut, &(self->timeStep) );
- Journal_Dump( self->forceOut, &(self->timeStep) );
-
KeyCall( self, self->loopElementsMomentumK, EntryPoint_VoidPtr_CallCast* )( KeyHandle( self, self->loopElementsMomentumK ), self );
+ _Snac_Context_WriteOutput( self );
+
+#if 0
/* Synchronise again, but really need to consolidate with the other one... */
if( self->forceCalcType == Snac_Force_Complete &&
self->mesh->layout->decomp->procsInUse > 1 )
{
Mesh_Sync( self->mesh );
}
-
- Journal_Dump( self->phaseIndexOut, &(self->timeStep) );
+#endif
}
@@ -1083,7 +1059,6 @@
void _Snac_Context_LoopNodes( void* context ) {
Snac_Context* self = (Snac_Context*)context;
Node_LocalIndex node_lI;
- Snac_Parallel* parallel = self->parallel;
if( self->rank == 0 ) Journal_DPrintf( self->debug, "In: %s\n", __func__ );
if( self->rank == 0 ) Journal_Printf(
@@ -1179,7 +1154,7 @@
MPI_Allreduce( &self->minLengthScale, &tmp, 1, MPI_DOUBLE, MPI_MIN, self->communicator );
self->minLengthScale = tmp;
- if(self->timeStep==self->restartStep && self->restartStep==0) {
+ if(self->timeStep==self->restartTimestep && self->restartTimestep==0) {
FILE* fp;
char fname[PATH_MAX];
self->initMinLengthScale = self->minLengthScale;
@@ -1202,29 +1177,51 @@
Journal_Printf( self->debug, "new self->dt: %g L/Vmax = %g\n", self->dt,self->minLengthScale/vmax );
Journal_Printf( self->debug, "new speed of sound: %g\n", self->speedOfSound );
}
+
+}
- fflush( self->timeStepInfo );
- Stream_Flush( self->strainRateOut );
- Stream_Flush( self->stressOut );
- Stream_Flush( self->hydroPressureOut );
- Stream_Flush( self->phaseIndexOut );
- Stream_Flush( self->coordOut );
- Stream_Flush( self->velOut );
- Stream_Flush( self->forceOut );
- fflush( self->stressTensorOut );
- fflush( NULL );
+void _Snac_Context_WriteLoopInfo( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
+
+ if( isTimeToDump( self ) ) {
+ fprintf(stderr,"writing timeStep.* at %d\n",self->timeStep );
+ _Snac_Context_DumpLoopInfo( self );
+ }
+
+/* if( isTimeToCheckpoint( self ) ) */
+/* _Snac_Context_CheckpointLoopInfo( self ); */
+
}
-void DumpLoopInfo( void* context ) {
+void _Snac_Context_DumpLoopInfo( void* context ) {
Snac_Context* self = (Snac_Context*)context;
fprintf( self->timeStepInfo, "%16u %16g %16g\n", self->timeStep, self->currentTime, self->dt );
+ fflush( self->timeStepInfo );
}
-void _Snac_Context_InitDump( Snac_Context* self ) {
- Snac_Mesh* tmpMesh;
+
+void _Snac_Context_CheckpointLoopInfo( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
+
+ fprintf( self->checkpointTimeStepInfo, "%16u %16g %16g\n", self->timeStep, self->currentTime, self->dt );
+ fflush( self->checkpointTimeStepInfo );
+}
+
+
+void _Snac_Context_InitOutput( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
+
+ _Snac_Context_InitDump( self );
+ _Snac_Context_InitCheckpoint( self );
+
+}
+
+
+void _Snac_Context_InitDump( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
char tmpBuf[200];
/* Create the strain rate dumping stream */
@@ -1304,81 +1301,168 @@
}
}
-void _Snac_Context_AdjustDump( Snac_Context* self ) {
- Snac_Mesh* tmpMesh;
+
+void _Snac_Context_InitCheckpoint( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
char tmpBuf[200];
- /* Adjust the strain rate dumping stream */
- sprintf( tmpBuf, "%s/strainRate.%u", self->outputPath, self->rank );
+ /* Create the phase index dumping stream */
+ self->phaseIndexCheckpoint = Journal_Register( VariableDumpStream_Type, "PhaseIndexCP" );
+ sprintf( tmpBuf, "%s/phaseIndexCP.%u", self->outputPath, self->rank );
VariableDumpStream_SetVariable(
- self->strainRateOut,
- Variable_Register_GetByName( self->variable_Register, "strainRate" ),
+ self->phaseIndexCheckpoint,
+ Variable_Register_GetByName( self->variable_Register, "elementMaterial" ),
self->mesh->elementLocalCount,
- self->dumpEvery,
+ self->checkpointEvery,
tmpBuf );
- /* Adjust the stress dumping stream */
- sprintf( tmpBuf, "%s/stress.%u", self->outputPath, self->rank );
+ /* Create the coords dumping stream */
+ self->coordCheckpoint = Journal_Register( VariableDumpStream_Type, "CoordCP" );
+ sprintf( tmpBuf, "%s/coordCP.%u", self->outputPath, self->rank );
VariableDumpStream_SetVariable(
- self->stressOut,
- Variable_Register_GetByName( self->variable_Register, "stress" ),
- self->mesh->elementLocalCount,
- self->dumpEvery,
- tmpBuf );
-
- /* Adjust the pressure dumping stream */
- sprintf( tmpBuf, "%s/hydroPressure.%u", self->outputPath, self->rank );
- VariableDumpStream_SetVariable(
- self->hydroPressureOut,
- Variable_Register_GetByName( self->variable_Register, "hydroPressure" ),
- self->mesh->elementLocalCount,
- self->dumpEvery,
- tmpBuf );
-
- /* Adjust the coords dumping stream */
- sprintf( tmpBuf, "%s/coord.%u", self->outputPath, self->rank );
- VariableDumpStream_SetVariable(
- self->coordOut,
+ self->coordCheckpoint,
Variable_Register_GetByName( self->variable_Register, "coord" ),
self->mesh->nodeLocalCount,
- self->dumpEvery,
+ self->checkpointEvery,
tmpBuf );
- /* Adjust the velocity dumping stream */
- sprintf( tmpBuf, "%s/vel.%u", self->outputPath, self->rank );
+ /* Create the velocity dumping stream */
+ self->velCheckpoint = Journal_Register( VariableDumpStream_Type, "VelocityCP" );
+ sprintf( tmpBuf, "%s/velCP.%u", self->outputPath, self->rank );
VariableDumpStream_SetVariable(
- self->velOut,
+ self->velCheckpoint,
Variable_Register_GetByName( self->variable_Register, "velocity" ),
self->mesh->nodeLocalCount,
- self->dumpEvery,
+ self->checkpointEvery,
tmpBuf );
- /* Adjust the force dumping stream */
- sprintf( tmpBuf, "%s/force.%u", self->outputPath, self->rank );
- VariableDumpStream_SetVariable(
- self->forceOut,
- Variable_Register_GetByName( self->variable_Register, "force" ),
- self->mesh->nodeLocalCount,
- self->dumpEvery,
- tmpBuf );
+ /* Create the stressTensor dump file */
+ sprintf( tmpBuf, "%s/stressTensorCP.%u", self->outputPath, self->rank );
+ if( (self->stressTensorCheckpoint = fopen( tmpBuf, "w+" )) == NULL ) {
+ assert( self->stressTensorCheckpoint /* failed to open file for writing */ );
+ abort();
+ }
}
-void _Snac_Context_DumpStressTensor( Snac_Context* self ) {
+void _Snac_Context_WriteOutput( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
+
+ if( isTimeToDump( self ) )
+ _Snac_Context_Dump( self );
+
+ if( isTimeToCheckpoint( self ) )
+ _Snac_Context_Checkpoint( self );
+}
+
+
+void _Snac_Context_Dump( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
- /* fprintf(stderr, "Dumping ? stress tensor: ts=%d df=%d\n",self->timeStep,self->dumpEvery); */
- if( self->timeStep ==0 || (self->timeStep-1) % self->dumpEvery == 0 ) {
- Element_LocalIndex element_lI;
+ Journal_Dump( self->strainRateOut, NULL );
+ Journal_Dump( self->stressOut, NULL );
+ Journal_Dump( self->hydroPressureOut, NULL );
+ _Snac_Context_DumpStressTensor( self );
+ Journal_Dump( self->coordOut, NULL );
+ Journal_Dump( self->velOut, NULL );
+ Journal_Dump( self->forceOut, NULL );
+ Journal_Dump( self->phaseIndexOut, NULL );
+
+ Stream_Flush( self->strainRateOut );
+ Stream_Flush( self->stressOut );
+ Stream_Flush( self->hydroPressureOut );
+ Stream_Flush( self->phaseIndexOut );
+ Stream_Flush( self->coordOut );
+ Stream_Flush( self->velOut );
+ Stream_Flush( self->forceOut );
+ fflush( self->stressTensorOut );
+ fflush( NULL );
+}
+
+
+void _Snac_Context_Checkpoint( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
+
+ _Snac_Context_CheckpointStressTensor( self );
+ Journal_Dump( self->coordCheckpoint, NULL );
+ Journal_Dump( self->velCheckpoint, NULL );
+ Journal_Dump( self->phaseIndexCheckpoint, NULL );
+
+ fflush( self->stressTensorCheckpoint );
+ Stream_Flush( self->coordCheckpoint );
+ Stream_Flush( self->velCheckpoint );
+ Stream_Flush( self->phaseIndexCheckpoint );
+}
+
+
+void _Snac_Context_DumpStressTensor( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
+
+ Element_LocalIndex element_lI;
+
+ for( element_lI = 0; element_lI < self->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Element* element = Snac_Element_At( self, element_lI );
+ Tetrahedra_Index tetra_I;
+ float stressVector[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+ float totalVolume = 0.0;
- for( element_lI = 0; element_lI < self->mesh->elementLocalCount; element_lI++ ) {
- Snac_Element* element = Snac_Element_At( self, element_lI );
- /* Take average of tetra viscosity for the element */
- Tetrahedra_Index tetra_I;
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- float stressVector[6] = { element->tetra[tetra_I].stress[0][0], element->tetra[tetra_I].stress[1][1], element->tetra[tetra_I].stress[2][2],
- element->tetra[tetra_I].stress[0][1], element->tetra[tetra_I].stress[0][2], element->tetra[tetra_I].stress[1][2]};
- fwrite( &stressVector, sizeof(float), 6, self->stressTensorOut );
- }
+ /* Take average of tetra stress for the element */
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ )
+ totalVolume += element->tetra[tetra_I].volume;
+
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ stressVector[0] += element->tetra[tetra_I].stress[0][0]*element->tetra[tetra_I].volume/totalVolume;
+ stressVector[1] += element->tetra[tetra_I].stress[1][1]*element->tetra[tetra_I].volume/totalVolume;
+ stressVector[2] += element->tetra[tetra_I].stress[2][2]*element->tetra[tetra_I].volume/totalVolume;
+ stressVector[3] += element->tetra[tetra_I].stress[0][1]*element->tetra[tetra_I].volume/totalVolume;
+ stressVector[4] += element->tetra[tetra_I].stress[0][2]*element->tetra[tetra_I].volume/totalVolume;
+ stressVector[5] += element->tetra[tetra_I].stress[1][2]*element->tetra[tetra_I].volume/totalVolume;
}
- }
+ fwrite( &stressVector, sizeof(float), 6, self->stressTensorOut );
+ }
}
+
+
+void _Snac_Context_CheckpointStressTensor( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
+
+ Element_LocalIndex element_lI;
+
+ for( element_lI = 0; element_lI < self->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Element* element = Snac_Element_At( self, element_lI );
+ /* Take average of tetra viscosity for the element */
+ Tetrahedra_Index tetra_I;
+
+ /* Write the stress vector for all the tets. This is for restarting. */
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ float stressVector[6] = { element->tetra[tetra_I].stress[0][0], element->tetra[tetra_I].stress[1][1], element->tetra[tetra_I].stress[2][2],
+ element->tetra[tetra_I].stress[0][1], element->tetra[tetra_I].stress[0][2], element->tetra[tetra_I].stress[1][2]};
+ fwrite( &stressVector, sizeof(float), 6, self->stressTensorCheckpoint );
+ }
+ }
+}
+
+
+Bool isTimeToDump( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
+
+ if( (self->timeStep==1) || ((self->timeStep>1)&&(self->timeStep%self->dumpEvery==0)) )
+ return True;
+ else
+ return False;
+
+}
+
+
+Bool isTimeToCheckpoint( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
+
+ if( (self->timeStep==self->maxTimeSteps) )
+ return True;
+ else {
+ if( self->checkpointEvery == 0 )
+ return False;
+ else if( (self->timeStep%self->checkpointEvery==0) )
+ return True;
+ }
+}
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Context.h 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Context.h 2009-04-10 20:43:44 UTC (rev 14656)
@@ -73,7 +73,6 @@
/* Snac_Context info */ \
Bool spherical; \
Bool computeThermalStress; \
- Index restartStep; \
\
\
double topo_kappa; \
@@ -100,14 +99,21 @@
/* TODO... we want journal or the like to look after this in the end */ \
FILE* simInfo; \
FILE* timeStepInfo; \
- FILE* phaseIndexOut; \
+ FILE* stressTensorOut; \
+ Stream* phaseIndexOut; \
Stream* strainRateOut; \
Stream* hydroPressureOut; \
Stream* stressOut; \
- FILE* stressTensorOut; \
Stream* coordOut; \
Stream* velOut; \
Stream* forceOut; \
+ \
+ FILE* checkpointTimeStepInfo; \
+ FILE* stressTensorCheckpoint; \
+ Stream* phaseIndexCheckpoint; \
+ Stream* coordCheckpoint; \
+ Stream* velCheckpoint; \
+ \
Stream* snacInfo; \
Stream* snacDebug; \
Stream* snacVerbose; \
@@ -223,16 +229,19 @@
void _Snac_Context_Sync( void* context );
/* Some output dumping helpers */
- void DumpLoopInfo( void* context );
- void DumpStrainRateAndStress( Snac_Context* self, Element_LocalIndex element_lI );
- void DumpCoord( Snac_Context* self, Node_LocalIndex node_lI );
- void DumpVelocity( Snac_Context* self, Node_LocalIndex node_lI );
- void DumpForces( Snac_Context* self, Node_LocalIndex node_lI );
- void DumpForceCheckSum( Snac_Context* self );
+ Bool isTimeToDump( void* context );
+ Bool isTimeToCheckpoint( void* context );
- void _Snac_Context_InitDump( Snac_Context* self );
- void _Snac_Context_AdjustDump( Snac_Context* self );
- void _Snac_Context_DumpStressTensor( Snac_Context* self );
- void _Snac_Context_DumpPhaseIndex( Snac_Context* self );
+ void _Snac_Context_WriteLoopInfo( void* context );
+ void _Snac_Context_DumpLoopInfo( void* context );
+ void _Snac_Context_CheckpointLoopInfo( void* context );
+ void _Snac_Context_InitOutput( void* context );
+ void _Snac_Context_InitDump( void* context );
+ void _Snac_Context_InitCheckpoint( void* context );
+ void _Snac_Context_WriteOutput( void* context );
+ void _Snac_Context_Dump( void* context );
+ void _Snac_Context_Checkpoint( void* context );
+ void _Snac_Context_DumpStressTensor( void* context );
+ void _Snac_Context_CheckpointStressTensor( void* context );
#endif /* __Snac_Context_h__ */
Modified: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Track.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -106,7 +106,7 @@
}
}
/* if( restart ) { */
-/* fprintf(stderr, "Restarting: in Track.c with ts=%d, restart ts=%d\n", context->timeStep, context->restartStep); */
+/* fprintf(stderr, "Restarting: in Track.c with ts=%d, restart ts=%d\n", context->timeStep, context->restartTimestep); */
/* return; */
/* } */
@@ -115,9 +115,9 @@
/*
* 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 || (restart && (context->timeStep-context->restartStep)==1)) {
+ if(context->timeStep==1 || (restart && (context->timeStep-context->restartTimestep)==1)) {
#ifdef DEBUG4
- fprintf(stderr,"Tracking: creating surface grid record: %d\n",context->timeStep-context->restartStep);
+ fprintf(stderr,"Tracking: creating surface grid record: %d\n",context->timeStep-context->restartTimestep);
#endif
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)));
@@ -170,7 +170,7 @@
*/
node_yVelocity = node_yElevation-*tmp_yGridOldPtr;
node_yAcceln = node_yVelocity-(*tmp_yGridOldPtr-*tmp_yGridOlderPtr);
- if(context->timeStep>=3 || (restart && (context->timeStep-context->restartStep)>=3)){
+ if(context->timeStep>=3 || (restart && (context->timeStep-context->restartTimestep)>=3)){
if(fabs(node_yVelocity)>max_yVelocity)
max_yVelocity = fabs(node_yVelocity);
if(fabs(node_yAcceln)>max_yAcceln)
@@ -211,7 +211,7 @@
/* if(unit_yAcceln==0.0 && max_yAcceln>0.0) */
/* unit_yAcceln = max_yAcceln; */
if( !contextExt->startedTrackingFlag && max_yVelocity>=startThreshold
- && (context->timeStep>=4 || (restart && (context->timeStep-context->restartStep)>=4)) )
+ && (context->timeStep>=4 || (restart && (context->timeStep-context->restartTimestep)>=4)) )
contextExt->startedTrackingFlag=TRUE;
#ifdef DEBUG
@@ -227,7 +227,7 @@
if( !contextExt->elasticStabilizedFlag
&& CheckStabilizingFn(max_yVelocity, max_yAcceln, stopThreshold, fallingFlag)==TRUE
&& (context->maxTimeSteps!=context->timeStep
- || (restart && (context->maxTimeSteps!=(context->timeStep-context->restartStep)) )) ) {
+ || (restart && (context->maxTimeSteps!=(context->timeStep-context->restartTimestep)) )) ) {
/*
* Stabilizing on this thread
*/
@@ -268,7 +268,7 @@
* In addition, force a dump of this model state by changing dump freq to 1.
*/
dumpEvery=1;
- maxTimeSteps=(!restart ? context->timeStep+1 : context->timeStep-context->restartStep+1);
+ maxTimeSteps=(!restart ? context->timeStep+1 : context->timeStep-context->restartTimestep+1);
doneTrackingFlag=TRUE;
#ifdef DEBUG4
fprintf(stderr,"r=%d, ts=%d/%d: Terminating simulation at step=%d\n",context->rank, context->timeStep,
@@ -281,7 +281,7 @@
dumpEvery=contextExt->plasticDeformationDumpFreq;
#ifdef DEBUG4
fprintf(stderr,"r=%d, ts=%d/%d: Shifting to elastoplastic simulation at step=%d for %d more steps\n",context->rank, context->timeStep,
- context->maxTimeSteps, context->timeStep, maxTimeSteps-(context->timeStep-context->restartStep));
+ context->maxTimeSteps, context->timeStep, maxTimeSteps-(context->timeStep-context->restartTimestep));
#endif
}
}
Modified: long/3D/SNAC/trunk/Snac/plugins/maxwell/ConstructExtensions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/maxwell/ConstructExtensions.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/maxwell/ConstructExtensions.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -75,9 +75,15 @@
context->condFunc_Register,
ConditionFunction_New( _SnacVelocity_VariableCondition, "variablevelBC" ) );
- /* Prepare the dump file */
+ /* Prepare the dump and checkpoint file */
sprintf( tmpBuf, "%s/viscosity.%u", context->outputPath, context->rank );
if( (contextExt->viscosityOut = fopen( tmpBuf, "w+" )) == NULL ) {
assert( contextExt->viscosityOut /* failed to open file for writing */ );
+ abort();
}
+ sprintf( tmpBuf, "%s/viscosityCP.%u", context->outputPath, context->rank );
+ if( (contextExt->viscosityCheckpoint = fopen( tmpBuf, "w+" )) == NULL ) {
+ assert( contextExt->viscosityCheckpoint /* failed to open file for writing */ );
+ abort();
+ }
}
Modified: long/3D/SNAC/trunk/Snac/plugins/maxwell/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/maxwell/Context.h 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/maxwell/Context.h 2009-04-10 20:43:44 UTC (rev 14656)
@@ -45,6 +45,7 @@
/* If we need to store anything on the context for viscosity... do it here */
/* Dumping */
FILE* viscosityOut;
+ FILE* viscosityCheckpoint;
};
/* Print the contents of the context extension */
Modified: long/3D/SNAC/trunk/Snac/plugins/maxwell/Output.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/maxwell/Output.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/maxwell/Output.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -42,38 +42,78 @@
#include <math.h>
#include <assert.h>
-void _SnacMaxwell_DumpViscosity( void* _context ) {
+void _SnacMaxwell_WriteViscosity( void* _context ) {
Snac_Context* context = (Snac_Context*) _context;
- SnacMaxwell_Context* contextExt = ExtensionManager_Get(
- context->extensionMgr,
- context,
- SnacMaxwell_ContextHandle );
- if( context->timeStep ==0 || (context->timeStep-1) % context->dumpEvery == 0 ) {
- Element_LocalIndex element_lI;
+ if( isTimeToDump( context ) )
+ _SnacMaxwell_DumpViscosity( context );
+ if( isTimeToCheckpoint( context ) )
+ _SnacMaxwell_CheckpointViscosity( context );
+}
- #if DEBUG
- printf( "In %s()\n", __func__ );
- #endif
- for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
- Snac_Element* element = Snac_Element_At( context, element_lI );
- SnacMaxwell_Element* elementExt = ExtensionManager_Get(
- context->mesh->elementExtensionMgr,
- element,
- SnacMaxwell_ElementHandle );
- /* Take average of tetra viscosity for the element */
- Tetrahedra_Index tetra_I;
- float viscosity = 0.0f;
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ )
- viscosity += elementExt->viscosity[tetra_I]/Tetrahedra_Count;
- assert(viscosity >= 0.0);
+void _SnacMaxwell_DumpViscosity( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ SnacMaxwell_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacMaxwell_ContextHandle );
+ Element_LocalIndex element_lI;
- if(viscosity > 0.0)
- viscosity = log10(viscosity);
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Element* element = Snac_Element_At( context, element_lI );
+ SnacMaxwell_Element* elementExt = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ element,
+ SnacMaxwell_ElementHandle );
+ /* Take average of tetra viscosity for the element */
+ Tetrahedra_Index tetra_I;
+ float viscosity = 0.0f;
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ )
+ viscosity += elementExt->viscosity[tetra_I]/Tetrahedra_Count;
+ assert(viscosity >= 0.0);
+
+ if(viscosity > 0.0)
+ viscosity = log10(viscosity);
+
+ fwrite( &viscosity, sizeof(float), 1, contextExt->viscosityOut );
+ }
+ fflush( contextExt->viscosityOut );
+}
- fwrite( &viscosity, sizeof(float), 1, contextExt->viscosityOut );
- }
- fflush( contextExt->viscosityOut );
+void _SnacMaxwell_CheckpointViscosity( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ SnacMaxwell_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacMaxwell_ContextHandle );
+ Element_LocalIndex element_lI;
+
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Element* element = Snac_Element_At( context, element_lI );
+ SnacMaxwell_Element* elementExt = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ element,
+ SnacMaxwell_ElementHandle );
+ /* Take average of tetra viscosity for the element */
+ Tetrahedra_Index tetra_I;
+ float viscosity = 0.0f;
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ )
+ viscosity += elementExt->viscosity[tetra_I]/Tetrahedra_Count;
+ assert(viscosity >= 0.0);
+
+ if(viscosity > 0.0)
+ viscosity = log10(viscosity);
+
+ fwrite( &viscosity, sizeof(float), 1, contextExt->viscosityCheckpoint );
}
+ fflush( contextExt->viscosityCheckpoint );
}
Modified: long/3D/SNAC/trunk/Snac/plugins/maxwell/Output.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/maxwell/Output.h 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/maxwell/Output.h 2009-04-10 20:43:44 UTC (rev 14656)
@@ -32,6 +32,8 @@
#ifndef _SnacMaxwell_Output_
#define _SnacMaxwell_Output_
+ void _SnacMaxwell_WriteViscosity( void* _context );
void _SnacMaxwell_DumpViscosity( void* _context );
+ void _SnacMaxwell_CheckpointViscosity( void* _context );
#endif
Modified: long/3D/SNAC/trunk/Snac/plugins/maxwell/Register.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/maxwell/Register.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/maxwell/Register.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -125,13 +125,13 @@
SnacMaxwell_Type );
EntryPoint_Prepend( /* Dump the initial viscosity */
Context_GetEntryPoint( context, AbstractContext_EP_Execute ),
- "SnacMaxwell_Dump",
- _SnacMaxwell_DumpViscosity,
+ "SnacMaxwell_Write",
+ _SnacMaxwell_WriteViscosity,
SnacMaxwell_Type );
EntryPoint_Append( /* and dump each loop */
Context_GetEntryPoint( context, Snac_EP_CalcStresses ),
- "SnacMaxwell_Dump",
- _SnacMaxwell_DumpViscosity,
+ "SnacMaxwell_Write",
+ _SnacMaxwell_WriteViscosity,
SnacMaxwell_Type );
/* Construct. */
Modified: long/3D/SNAC/trunk/Snac/plugins/plastic/Constitutive.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/plastic/Constitutive.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/plastic/Constitutive.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -193,7 +193,7 @@
if( h < 0.0f ) {
/* !shear failure */
- alam = fs / ( a1 - a2 * anpsi + a1 * anphi * anpsi - a2 * anphi + hardening );
+ alam = fs / ( a1 - a2 * anpsi + a1 * anphi * anpsi - a2 * anphi + 2.0*sqrt(anphi)*hardening );
s[0] -= alam * ( a1 - a2 * anpsi );
s[1] -= alam * a2 * ( 1.0f - anpsi );
s[2] -= alam * ( a2 - a1 * anpsi );
Modified: long/3D/SNAC/trunk/Snac/plugins/plastic/ConstructExtensions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/plastic/ConstructExtensions.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/plastic/ConstructExtensions.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -43,6 +43,30 @@
#define PATH_MAX 1024
#endif
+void _SnacVelocity_VariableCondition( Index node_dI, Variable_Index var_I, void* _context, void* result ){
+ Snac_Context* context = (Snac_Context*)_context;
+ Mesh* mesh = context->mesh;
+ MeshLayout* layout = (MeshLayout*)mesh->layout;
+ HexaMD* decomp = (HexaMD*)layout->decomp;
+
+ double* velComponent = (double*)result;
+ IJK ijk;
+ Node_GlobalIndex node_gI = _MeshDecomp_Node_LocalToGlobal1D( decomp, node_dI );
+ Index oneThird = (unsigned int)((11.0/30.0)*decomp->nodeGlobal3DCounts[0]);
+ Index twoThirds = (unsigned int)((19.0/30.0)*decomp->nodeGlobal3DCounts[0]);
+
+ const double vmag = 1.4e-06;
+
+ RegularMeshUtils_Node_1DTo3D( decomp, node_gI, &ijk[0], &ijk[1], &ijk[2] );
+
+ if( ijk[0] < oneThird )
+ (*velComponent) = -1.0*vmag;
+ else if( ijk[0] >= oneThird && ijk[0] <= twoThirds )
+ (*velComponent) = -1.0*vmag + 2.0*vmag*((double)(ijk[0])-(double)oneThird)/((double)twoThirds-(double)oneThird);
+ else
+ (*velComponent) = vmag;
+}
+
void _SnacPlastic_ConstructExtensions( void* _context, void* data ) {
Snac_Context* context = (Snac_Context*)_context;
SnacPlastic_Context* contextExt = ExtensionManager_Get(
@@ -55,9 +79,19 @@
printf( "In %s()\n", __func__ );
#endif
- /* Prepare the dump file */
+ ConditionFunction_Register_Add(
+ context->condFunc_Register,
+ ConditionFunction_New( _SnacVelocity_VariableCondition, "variableVelBC" ) );
+
+ /* Prepare the dump and checkpoint file */
sprintf( tmpBuf, "%s/plStrain.%u", context->outputPath, context->rank );
if( (contextExt->plStrainOut = fopen( tmpBuf, "w+" )) == NULL ) {
assert( contextExt->plStrainOut /* failed to open file for writing */ );
+ abort();
}
+ sprintf( tmpBuf, "%s/plStrainCP.%u", context->outputPath, context->rank );
+ if( (contextExt->plStrainCheckpoint = fopen( tmpBuf, "w+" )) == NULL ) {
+ assert( contextExt->plStrainCheckpoint /* failed to open file for writing */ );
+ abort();
+ }
}
Modified: long/3D/SNAC/trunk/Snac/plugins/plastic/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/plastic/Context.h 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/plastic/Context.h 2009-04-10 20:43:44 UTC (rev 14656)
@@ -44,6 +44,7 @@
struct _SnacPlastic_Context {
/* Dumping */
FILE* plStrainOut;
+ FILE* plStrainCheckpoint;
};
/* Print the contents of the context extension */
Modified: long/3D/SNAC/trunk/Snac/plugins/plastic/Output.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/plastic/Output.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/plastic/Output.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -40,31 +40,69 @@
#include "Register.h"
#include <stdio.h>
+
+
+void _SnacPlastic_WritePlasticStrain( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+
+ if( isTimeToDump( context ) )
+ _SnacPlastic_DumpPlasticStrain( context );
+
+ if( isTimeToCheckpoint( context ) )
+ _SnacPlastic_CheckpointPlasticStrain( context );
+
+}
+
+
void _SnacPlastic_DumpPlasticStrain( void* _context ) {
Snac_Context* context = (Snac_Context*) _context;
- SnacPlastic_Context* contextExt = ExtensionManager_Get(
- context->extensionMgr,
- context,
- SnacPlastic_ContextHandle );
+ SnacPlastic_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacPlastic_ContextHandle );
+ Element_LocalIndex element_lI;
+
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Element* element = Snac_Element_At( context, element_lI );
+ SnacPlastic_Element* elementExt = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ element,
+ SnacPlastic_ElementHandle );
+ float plasticStrain = elementExt->aps;
+ /* Take average of tetra plastic strain for the element */
+ fwrite( &plasticStrain, sizeof(float), 1, contextExt->plStrainOut );
+ }
+ fflush( contextExt->plStrainOut );
+}
- if( context->timeStep ==0 || (context->timeStep-1) % context->dumpEvery == 0 ) {
- Element_LocalIndex element_lI;
+void _SnacPlastic_CheckpointPlasticStrain( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ SnacPlastic_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacPlastic_ContextHandle );
+ Element_LocalIndex element_lI;
- #if DEBUG
- printf( "In %s()\n", __func__ );
- #endif
-
- for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
- Snac_Element* element = Snac_Element_At( context, element_lI );
- SnacPlastic_Element* elementExt = ExtensionManager_Get(
- context->mesh->elementExtensionMgr,
- element,
- SnacPlastic_ElementHandle );
- float plasticStrain = elementExt->aps;
- /* Take average of tetra plastic strain for the element */
- fwrite( &plasticStrain, sizeof(float), 1, contextExt->plStrainOut );
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Element* element = Snac_Element_At( context, element_lI );
+ SnacPlastic_Element* plasticElement = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ element,
+ SnacPlastic_ElementHandle );
+ Tetrahedra_Index tetra_I;
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ float plasticStrain = plasticElement->plasticStrain[tetra_I];
+ fwrite( &plasticStrain, sizeof(float), 1, contextExt->plStrainCheckpoint );
}
- fflush( contextExt->plStrainOut );
}
+ fflush( contextExt->plStrainCheckpoint );
}
Modified: long/3D/SNAC/trunk/Snac/plugins/plastic/Output.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/plastic/Output.h 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/plastic/Output.h 2009-04-10 20:43:44 UTC (rev 14656)
@@ -32,6 +32,8 @@
#ifndef _SnacPlastic_Output_
#define _SnacPlastic_Output_
+ void _SnacPlastic_WritePlasticStrain( void* _context );
void _SnacPlastic_DumpPlasticStrain( void* _context );
+ void _SnacPlastic_CheckpointPlasticStrain( void* _context );
#endif
Modified: long/3D/SNAC/trunk/Snac/plugins/plastic/Register.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/plastic/Register.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/plastic/Register.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -109,13 +109,13 @@
SnacPlastic_Type );
EntryPoint_Prepend( /* Dump the initial plastic strain */
Context_GetEntryPoint( context, AbstractContext_EP_Execute ),
- "SnacPlastic_Dump",
- _SnacPlastic_DumpPlasticStrain,
+ "SnacPlastic_Write",
+ _SnacPlastic_WritePlasticStrain,
SnacPlastic_Type );
EntryPoint_Append( /* and dump each loop */
Context_GetEntryPoint( context, Snac_EP_CalcStresses ),
- "SnacPlastic_Dump",
- _SnacPlastic_DumpPlasticStrain,
+ "SnacPlastic_Write",
+ _SnacPlastic_WritePlasticStrain,
SnacPlastic_Type );
/* Add extensions to the interpolate element entry point, but it will only exist if the remesher is loaded. */
Modified: long/3D/SNAC/trunk/Snac/plugins/plastic/Remesh.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/plastic/Remesh.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/plastic/Remesh.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -65,4 +65,9 @@
SnacPlastic_ElementHandle );
elementExt->plasticStrain[dstTetInd] = fromElementExt->plasticStrain[srcTetInd];
+
+ if( (context->rank==1 && (dstElementInd <=7 && dstElementInd >= 6) ) ) //(context->rank==2 && dstElementInd <=2 )
+ fprintf(stderr, "element_lI: (%u %u %e), fromElement_lI: (%u %u %e)\n",
+ dstElementInd, dstTetInd, elementExt->plasticStrain[dstTetInd],
+ srcElementInd, srcTetInd, fromElementExt->plasticStrain[srcTetInd] );
}
Modified: long/3D/SNAC/trunk/Snac/plugins/plastic/makefile
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/plastic/makefile 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/plastic/makefile 2009-04-10 20:43:44 UTC (rev 14656)
@@ -48,6 +48,6 @@
EXTERNAL_LIBS = -L${STGERMAIN_LIBDIR} -lSnac -lStGermain
EXTERNAL_INCLUDES = -I${STGERMAIN_INCDIR}/StGermain -I${STGERMAIN_INCDIR}
-packages = MPI XML MATH
+packages = MPI XML MATH GSL
include ${PROJ_ROOT}/Makefile.vmake
Modified: long/3D/SNAC/trunk/Snac/plugins/remesher/InitialConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/remesher/InitialConditions.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/remesher/InitialConditions.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -77,12 +77,12 @@
/*
** Firstly, backup all coords into the 'newNodeCoords' coord block.
*/
- if( context->restartStep > 0 ) {
+ if( context->restartTimestep > 0 ) {
FILE *fp;
char path[PATH_MAX];
double x,y,z;
- sprintf(path, "%s/snac.initcoord.%d.%06d.restart",context->outputPath,context->rank,context->restartStep);
+ sprintf(path, "%s/snac.initcoord.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
Journal_Firewall( ( (fp=fopen(path,"r")) != NULL), context->snacError, "Can't find %s", path );
for( lNode_i = 0; lNode_i < mesh->nodeLocalCount; lNode_i++ ) {
fscanf( fp, "%le %le %le", &x,&y,&z);
Modified: long/3D/SNAC/trunk/Snac/plugins/restarter/InitialConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/restarter/InitialConditions.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/restarter/InitialConditions.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -66,11 +66,11 @@
Journal_Printf( context->snacInfo, "In: %s\n", __func__ );
- sprintf(path, "%s/snac.coord.%d.%06d.restart",context->outputPath,context->rank,context->restartStep);
+ sprintf(path, "%s/snac.coord.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
#ifdef DEBUG
- fprintf(stderr,"Restarter: loading mesh coordinates from %s for r=%d, ts=%d \n",path,context->rank,context->restartStep);
+ fprintf(stderr,"Restarter: loading mesh coordinates from %s for r=%d, ts=%d \n",path,context->rank,context->restartTimestep);
#endif
- Journal_Firewall( (fp = fopen(path,"r")) != NULL, "Can't find %s - is the parameter \"restartStep\" set correctly in the input xml?\n", path);
+ Journal_Firewall( (fp = fopen(path,"r")) != NULL, "Can't find %s - is the parameter \"restartTimestep\" set correctly in the input xml?\n", path);
/* read in restart file to construct the initial mesh */
for( node_lI = 0; node_lI < context->mesh->nodeLocalCount; node_lI++ ) {
@@ -96,11 +96,11 @@
Journal_Printf( context->snacInfo, "In: %s\n", __func__ );
- sprintf(path, "%s/snac.vel.%d.%06d.restart",context->outputPath,context->rank,context->restartStep);
+ sprintf(path, "%s/snac.vel.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
#ifdef DEBUG
- fprintf(stderr,"Restarter: loading mesh velocities from %s for r=%d, ts=%d \n",path,context->rank,context->restartStep);
+ fprintf(stderr,"Restarter: loading mesh velocities from %s for r=%d, ts=%d \n",path,context->rank,context->restartTimestep);
#endif
- Journal_Firewall( (fp = fopen(path,"r")) != NULL, "Can't find %s - is the parameter \"restartStep\" set correctly in the input xml?\n", path);
+ Journal_Firewall( (fp = fopen(path,"r")) != NULL, "Can't find %s - is the parameter \"restartTimestep\" set correctly in the input xml?\n", path);
/* read in restart file to construct the initial mesh */
for( node_lI = 0; node_lI < context->mesh->nodeLocalCount; node_lI++ ) {
@@ -122,11 +122,11 @@
char path[PATH_MAX];
Journal_Printf( context->snacInfo, "In: %s\n", __func__ );
- sprintf(path, "%s/snac.stressTensor.%d.%06d.restart",context->outputPath,context->rank,context->restartStep);
- Journal_Firewall( (fp = fopen(path,"r")) != NULL, "Can't find %s - is the parameter \"restartStep\" set correctly in the input xml?\n", path );
+ sprintf(path, "%s/snac.stressTensor.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
+ Journal_Firewall( (fp = fopen(path,"r")) != NULL, "Can't find %s - is the parameter \"restartTimestep\" set correctly in the input xml?\n", path );
#ifdef DEBUG
- fprintf(stderr,"Restarter: loading stress tensors from %s for r=%d, ts=%d \n",path,context->rank,context->restartStep);
+ fprintf(stderr,"Restarter: loading stress tensors from %s for r=%d, ts=%d \n",path,context->rank,context->restartTimestep);
#endif
/* read in restart file to assign initial stress field. */
@@ -147,9 +147,9 @@
for(j=0;j<3;j++)
element->tetra[tetra_I].stress[i][j] = S[i][j];
}
- element->hydroPressure = 0.0f;
+ element->hydroPressure = 0.0;
for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- element->hydroPressure += element->tetra[tetra_I].stress[0][0] / Tetrahedra_Count;
+ element->hydroPressure += ((element->tetra[tetra_I].stress[0][0]+element->tetra[tetra_I].stress[1][1]+element->tetra[tetra_I].stress[2][2])/3.0/Tetrahedra_Count);
sVolAvg +=
element->tetra[tetra_I].stress[1][1] * element->tetra[tetra_I].stress[2][2] +
element->tetra[tetra_I].stress[2][2] * element->tetra[tetra_I].stress[0][0] +
@@ -161,7 +161,7 @@
}
sVolAvg /= Tetrahedra_Count;
sOtherAvg /= Tetrahedra_Count;
- element->stress = 0.5f * sqrt( 0.5f * fabs( -1.0f * sVolAvg + sOtherAvg ) );
+ element->stress = 0.5 * sqrt( 0.5 * fabs( -1.0f * sVolAvg + sOtherAvg ) );
}
fclose( fp );
}
Modified: long/3D/SNAC/trunk/Snac/plugins/restarter_old/InitialConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/restarter_old/InitialConditions.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/restarter_old/InitialConditions.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -66,11 +66,11 @@
Journal_Printf( context->snacInfo, "In: %s\n", __func__ );
- sprintf(path, "%s/snac.coord.%d.%06d.restart",context->outputPath,context->rank,context->restartStep);
+ sprintf(path, "%s/snac.coord.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
#ifdef DEBUG
- fprintf(stderr,"RestartOld: loading mesh coordinates from %s for r=%d, ts=%d \n",path,context->rank,context->restartStep);
+ fprintf(stderr,"RestartOld: loading mesh coordinates from %s for r=%d, ts=%d \n",path,context->rank,context->restartTimestep);
#endif
- Journal_Firewall( (fp = fopen(path,"r")) != NULL, "Can't find %s - is the parameter \"restartStep\" set correctly in the input xml?\n", path);
+ Journal_Firewall( (fp = fopen(path,"r")) != NULL, "Can't find %s - is the parameter \"restartTimestep\" set correctly in the input xml?\n", path);
/* read in restart file to construct the initial mesh */
for( node_lI = 0; node_lI < context->mesh->nodeLocalCount; node_lI++ ) {
@@ -96,11 +96,11 @@
Journal_Printf( context->snacInfo, "In: %s\n", __func__ );
- sprintf(path, "%s/snac.vel.%d.%06d.restart",context->outputPath,context->rank,context->restartStep);
+ sprintf(path, "%s/snac.vel.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
#ifdef DEBUG
- fprintf(stderr,"RestartOld: loading mesh velocities from %s for r=%d, ts=%d \n",path,context->rank,context->restartStep);
+ fprintf(stderr,"RestartOld: loading mesh velocities from %s for r=%d, ts=%d \n",path,context->rank,context->restartTimestep);
#endif
- Journal_Firewall( (fp = fopen(path,"r")) != NULL, "Can't find %s - is the parameter \"restartStep\" set correctly in the input xml?\n", path);
+ Journal_Firewall( (fp = fopen(path,"r")) != NULL, "Can't find %s - is the parameter \"restartTimestep\" set correctly in the input xml?\n", path);
/* read in restart file to construct the initial mesh */
for( node_lI = 0; node_lI < context->mesh->nodeLocalCount; node_lI++ ) {
@@ -122,11 +122,11 @@
char path[PATH_MAX];
Journal_Printf( context->snacInfo, "In: %s\n", __func__ );
- sprintf(path, "%s/snac.stressTensor.%d.%06d.restart",context->outputPath,context->rank,context->restartStep);
- Journal_Firewall( (fp = fopen(path,"r")) != NULL, "Can't find %s - is the parameter \"restartStep\" set correctly in the input xml?\n", path );
+ sprintf(path, "%s/snac.stressTensor.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
+ Journal_Firewall( (fp = fopen(path,"r")) != NULL, "Can't find %s - is the parameter \"restartTimestep\" set correctly in the input xml?\n", path );
#ifdef DEBUG
- fprintf(stderr,"RestartOld: loading stress tensors from %s for r=%d, ts=%d \n",path,context->rank,context->restartStep);
+ fprintf(stderr,"RestartOld: loading stress tensors from %s for r=%d, ts=%d \n",path,context->rank,context->restartTimestep);
#endif
/* read in restart file to assign initial stress field. */
@@ -144,9 +144,9 @@
for(j=0;j<3;j++)
element->tetra[tetra_I].stress[i][j] = S[i][j];
}
- element->hydroPressure = 0.0f;
+ element->hydroPressure = 0.0;
for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- element->hydroPressure += element->tetra[tetra_I].stress[0][0] / Tetrahedra_Count;
+ element->hydroPressure += ((element->tetra[tetra_I].stress[0][0]+element->tetra[tetra_I].stress[1][1]+element->tetra[tetra_I].stress[2][2])/3.0/Tetrahedra_Count);
sVolAvg +=
element->tetra[tetra_I].stress[1][1] * element->tetra[tetra_I].stress[2][2] +
element->tetra[tetra_I].stress[2][2] * element->tetra[tetra_I].stress[0][0] +
@@ -158,7 +158,7 @@
}
sVolAvg /= Tetrahedra_Count;
sOtherAvg /= Tetrahedra_Count;
- element->stress = 0.5f * sqrt( 0.5f * fabs( -1.0f * sVolAvg + sOtherAvg ) );
+ element->stress = 0.5 * sqrt( 0.5 * fabs( -1.0f * sVolAvg + sOtherAvg ) );
}
fclose( fp );
}
Modified: long/3D/SNAC/trunk/Snac/plugins/temperature/ConstructExtensions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/temperature/ConstructExtensions.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/temperature/ConstructExtensions.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -111,9 +111,15 @@
temperatureBCsDict,
context->mesh );
- /* Prepare the dump file */
+ /* Prepare the dump and checkpoint file */
sprintf( tmpBuf, "%s/temperature.%u", context->outputPath, context->rank );
if( (contextExt->temperatureOut = fopen( tmpBuf, "w+" )) == NULL ) {
assert( contextExt->temperatureOut /* failed to open file for writing */ );
+ abort();
}
+ sprintf( tmpBuf, "%s/temperatureCP.%u", context->outputPath, context->rank );
+ if( (contextExt->temperatureCheckpoint = fopen( tmpBuf, "w+" )) == NULL ) {
+ assert( contextExt->temperatureCheckpoint /* failed to open file for writing */ );
+ abort();
+ }
}
Modified: long/3D/SNAC/trunk/Snac/plugins/temperature/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/temperature/Context.h 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/temperature/Context.h 2009-04-10 20:43:44 UTC (rev 14656)
@@ -51,6 +51,7 @@
/* Dumping */
FILE* temperatureOut;
+ FILE* temperatureCheckpoint;
};
/* Print the contents of the context extension */
Modified: long/3D/SNAC/trunk/Snac/plugins/temperature/Output.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/temperature/Output.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/temperature/Output.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -40,30 +40,63 @@
#include "Register.h"
#include <stdio.h>
+void _SnacTemperature_WriteTemp( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+
+ if( isTimeToDump( context ) )
+ _SnacTemperature_DumpTemp( context );
+
+ if( isTimeToCheckpoint( context ) )
+ _SnacTemperature_CheckpointTemp( context );
+
+}
+
+
void _SnacTemperature_DumpTemp( void* _context ) {
Snac_Context* context = (Snac_Context*) _context;
- SnacTemperature_Context* contextExt = ExtensionManager_Get(
- context->extensionMgr,
- context,
- SnacTemperature_ContextHandle );
+ SnacTemperature_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacTemperature_ContextHandle );
Node_LocalIndex node_lI;
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ for( node_lI = 0; node_lI < context->mesh->nodeLocalCount; node_lI++ ) {
+ Snac_Node* node = Snac_Node_At( context, node_lI );
+ SnacTemperature_Node* nodeExt = ExtensionManager_Get(
+ context->mesh->nodeExtensionMgr,
+ node,
+ SnacTemperature_NodeHandle );
+ float temperature = nodeExt->temperature;
+ fwrite( &temperature, sizeof(float), 1, contextExt->temperatureOut );
+ }
+ fflush( contextExt->temperatureOut );
+}
- if( context->timeStep == 0 || (context->timeStep - 1) % context->dumpEvery == 0 ) {
- #if DEBUG
- printf( "In %s()\n", __func__ );
- #endif
- for( node_lI = 0; node_lI < context->mesh->nodeLocalCount; node_lI++ ) {
- Snac_Node* node = Snac_Node_At( context, node_lI );
- SnacTemperature_Node* nodeExt = ExtensionManager_Get(
+void _SnacTemperature_CheckpointTemp( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ SnacTemperature_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacTemperature_ContextHandle );
+ Node_LocalIndex node_lI;
+
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ for( node_lI = 0; node_lI < context->mesh->nodeLocalCount; node_lI++ ) {
+ Snac_Node* node = Snac_Node_At( context, node_lI );
+ SnacTemperature_Node* nodeExt = ExtensionManager_Get(
context->mesh->nodeExtensionMgr,
node,
SnacTemperature_NodeHandle );
- float temperature = nodeExt->temperature;
- fwrite( &temperature, sizeof(float), 1, contextExt->temperatureOut );
- }
- fflush( contextExt->temperatureOut );
+ float temperature = nodeExt->temperature;
+ fwrite( &temperature, sizeof(float), 1, contextExt->temperatureCheckpoint );
}
+ fflush( contextExt->temperatureCheckpoint );
}
-
Modified: long/3D/SNAC/trunk/Snac/plugins/temperature/Output.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/temperature/Output.h 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/temperature/Output.h 2009-04-10 20:43:44 UTC (rev 14656)
@@ -32,6 +32,8 @@
#ifndef _SnacTemperature_Output_
#define _SnacTemperature_Output_
+ void _SnacTemperature_WriteTemp( void* _context );
void _SnacTemperature_DumpTemp( void* _context );
+ void _SnacTemperature_CheckpointTemp( void* _context );
#endif
Modified: long/3D/SNAC/trunk/Snac/plugins/temperature/Register.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/temperature/Register.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/temperature/Register.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -137,13 +137,13 @@
SnacTemperature_Type );
EntryPoint_Prepend( /* Dump the initial temperature */
Context_GetEntryPoint( context, AbstractContext_EP_Execute ),
- "SnacTemperature_Dump",
- _SnacTemperature_DumpTemp,
+ "SnacTemperature_Write",
+ _SnacTemperature_WriteTemp,
SnacTemperature_Type );
EntryPoint_Append( /* and dump each loop */
Context_GetEntryPoint( context, Snac_EP_LoopNodesEnergy ),
- "SnacTemperature_Dump",
- _SnacTemperature_DumpTemp,
+ "SnacTemperature_Write",
+ _SnacTemperature_WriteTemp,
SnacTemperature_Type );
EntryPoint_Append(
Context_GetEntryPoint( context, AbstractContext_EP_DestroyExtensions ),
Modified: long/3D/SNAC/trunk/Snac/plugins/temperature/VariableConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/temperature/VariableConditions.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/temperature/VariableConditions.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -150,7 +150,7 @@
RegularMeshUtils_Node_1DTo3D( decomp, node_gI, &ijk[0], &ijk[1], &ijk[2] );
/* for Cartesian case */
- rMin = R - 4.0e+04f;
+ rMin = R - 3.0e+03f;
rMax = R;
r = rMax + (*coord)[1];
/*ccccc*/
@@ -160,9 +160,12 @@
rMin /= R;
rMax /= R;
+#if 0
if( ijk[0] >= midI-5 && ijk[0] <= midI+5 ) {
age = 0.01f + 1.0f*abs(ijk[0]-midI)/5.0*5.01;
}
+#endif
+ age = 0.5 - 0.2*ijk[2]/(decomp->nodeGlobal3DCounts[2]-1);
temp = (rMax-r) * 0.5f / sqrt(age/scalet);
*temperature = rTemp * erf(temp);
@@ -177,34 +180,20 @@
SnacTemperature_ContextHandle );
Element_LocalIndex element_lI;
- int restart = 0;
Dictionary_Entry_Value* pluginsList;
Dictionary_Entry_Value* plugin;
- pluginsList = PluginsManager_GetPluginsList( context->dictionary );
- if (pluginsList) {
- plugin = Dictionary_Entry_Value_GetFirstElement(pluginsList);
- while ( plugin ) {
- if ( 0 == strcmp( Dictionary_Entry_Value_AsString( plugin ),
- "SnacRestart" ) ) {
- restart = 1;
- break;
- }
- plugin = plugin->next;
- }
- }
-
Journal_Printf( context->snacInfo, "In: %s\n", __func__ );
/* Temperature ICs are applied onto the "nodeICs" of Snac, and hence do not need to be repeated, but we must apply BCs */
/* In case of restarting, Temperature IC is still applied here */
- if( restart ) {
+ if( context->restartTimestep>0 ) {
FILE* fp;
Node_LocalIndex node_lI;
char path[PATH_MAX];
- sprintf(path, "%s/snac.temp.%d.%06d.restart",context->outputPath,context->rank,context->restartStep);
+ sprintf(path, "%s/snac.temp.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
Journal_Firewall( ( (fp = fopen(path,"r")) != NULL ), context->snacError, "Cannot find %s", path );
/* read in restart file to construct the initial temperature */
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Constitutive.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Constitutive.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Constitutive.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -324,7 +324,7 @@
if( h < 0.0f ) {
/* !shear failure */
- alam = fs / ( a1 - a2 * anpsi + a1 * anphi * anpsi - a2 * anphi + hardening );
+ alam = fs / ( a1 - a2 * anpsi + a1 * anphi * anpsi - a2 * anphi + 2.0*sqrt(anphi)*hardening );
s[0] -= alam * ( a1 - a2 * anpsi );
s[1] -= alam * a2 * ( 1.0f - anpsi );
s[2] -= alam * ( a2 - a1 * anpsi );
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic/ConstructExtensions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic/ConstructExtensions.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic/ConstructExtensions.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -108,14 +108,21 @@
sprintf( tmpBuf, "%s/plStrain.%u", context->outputPath, context->rank );
if( (contextExt->plStrainOut = fopen( tmpBuf, "w+" )) == NULL ) {
assert( contextExt->plStrainOut /* failed to open file for writing */ );
+ abort();
}
- /* for restarting, store each Tet's plastic strain */
- sprintf( tmpBuf, "%s/plStrainTensor.%u", context->outputPath, context->rank );
- if( (contextExt->plstrainTensorOut = fopen( tmpBuf, "w+" )) == NULL ) {
- assert( contextExt->plstrainTensorOut /* failed to open file for writing */ );
+ sprintf( tmpBuf, "%s/plStrainCP.%u", context->outputPath, context->rank );
+ if( (contextExt->plStrainCheckpoint = fopen( tmpBuf, "w+" )) == NULL ) {
+ assert( contextExt->plStrainCheckpoint /* failed to open file for writing */ );
+ abort();
}
sprintf( tmpBuf, "%s/viscosity.%u", context->outputPath, context->rank );
if( (contextExt->viscOut = fopen( tmpBuf, "w+" )) == NULL ) {
assert( contextExt->viscOut /* failed to open file for writing */ );
+ abort();
}
+ sprintf( tmpBuf, "%s/viscosityCP.%u", context->outputPath, context->rank );
+ if( (contextExt->viscCheckpoint = fopen( tmpBuf, "w+" )) == NULL ) {
+ assert( contextExt->viscCheckpoint /* failed to open file for writing */ );
+ abort();
+ }
}
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Context.h 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Context.h 2009-04-10 20:43:44 UTC (rev 14656)
@@ -44,8 +44,10 @@
struct _SnacViscoPlastic_Context {
/* Dumping */
FILE* plStrainOut;
+ FILE* viscOut;
+ FILE* plStrainCheckpoint;
+ FILE* viscCheckpoint;
FILE* plstrainTensorOut;
- FILE* viscOut;
};
/* Print the contents of the context extension */
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic/InitialConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic/InitialConditions.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic/InitialConditions.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -43,7 +43,7 @@
#include <math.h>
#define min(a,b) ((a)<(b) ? (a):(b))
#ifndef PATH_MAX
- #define PATH_MAX 1024
+#define PATH_MAX 1024
#endif
void SnacViscoPlastic_InitialConditions( void* _context, void* data ) {
@@ -53,102 +53,65 @@
double dt_maxwell = 0;
const double mfrac = 9.0e-01;
- int restart = 0;
- Dictionary_Entry_Value* pluginsList;
- Dictionary_Entry_Value* plugin;
- if( context->rank == 0 ) Journal_Printf( context->snacInfo, "In: %s\n", __func__ );
-// search for restart plugin
- pluginsList = PluginsManager_GetPluginsList( context->dictionary );
- if (pluginsList) {
- plugin = Dictionary_Entry_Value_GetFirstElement(pluginsList);
- while ( plugin ) {
- if ( 0 == strcmp( Dictionary_Entry_Value_AsString( plugin ),
- "SnacRestart" ) ) {
- restart = 1;
- break;
- }
- plugin = plugin->next;
- }
- }
-
-// loop over the phase using the dictionary
+ // loop over the phase using the dictionary
- Dictionary_Entry_Value* materialList = Dictionary_Get( context->dictionary, "materials" );
- int PhaseI = 0;
- if( materialList ) {
- Dictionary_Entry_Value* materialEntry = Dictionary_Entry_Value_GetFirstElement( materialList );
- /* loop around the phases to initialize rheology */
- while( materialEntry ) {
- context->materialProperty[PhaseI].rheology |= Snac_Material_ViscoPlastic;
- mu = context->materialProperty[PhaseI].mu;
- vis_min = context->materialProperty[PhaseI].vis_min;
- dt_maxwellI = mfrac*vis_min/mu;
+ Dictionary_Entry_Value* materialList = Dictionary_Get( context->dictionary, "materials" );
+ int PhaseI = 0;
+ if( materialList ) {
+ Dictionary_Entry_Value* materialEntry = Dictionary_Entry_Value_GetFirstElement( materialList );
+ /* loop around the phases to initialize rheology */
+ while( materialEntry ) {
+ context->materialProperty[PhaseI].rheology |= Snac_Material_ViscoPlastic;
+ mu = context->materialProperty[PhaseI].mu;
+ vis_min = context->materialProperty[PhaseI].vis_min;
+ dt_maxwellI = mfrac*vis_min/mu;
dt_maxwell = (PhaseI==0)?dt_maxwellI:min(dt_maxwellI,dt_maxwell);
- PhaseI++;
- materialEntry = materialEntry->next;
+ PhaseI++;
+ materialEntry = materialEntry->next;
- }
- }
- else {
- context->materialProperty[PhaseI].rheology |= Snac_Material_ViscoPlastic;
- mu = context->materialProperty[PhaseI].mu;
- vis_min = context->materialProperty[PhaseI].vis_min;
- dt_maxwell = mfrac*vis_min/mu;
- }
+ }
+ }
+ else {
+ context->materialProperty[PhaseI].rheology |= Snac_Material_ViscoPlastic;
+ mu = context->materialProperty[PhaseI].mu;
+ vis_min = context->materialProperty[PhaseI].vis_min;
+ dt_maxwell = mfrac*vis_min/mu;
+ }
if( context->dt > dt_maxwell) {
if(context->dtType == Snac_DtType_Constant)
fprintf(stderr,"dt(%e) should be smaller than dt_maxwell(%e) (mu=%e vis_min=%e mfrac=%e)\n",context->dt,dt_maxwell,mu,vis_min,mfrac);
else context->dt = dt_maxwell;
- }
- if( restart ) {
+ }
+ if( context->restartTimestep > 0 ) {
FILE* fp1;
FILE* fp2;
char path[PATH_MAX];
- sprintf(path, "%s/snac.plStrainTensor.%d.%06d.restart",context->outputPath,context->rank,context->restartStep);
+ sprintf(path, "%s/snac.plStrainTensor.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
Journal_Firewall( (fp1 = fopen(path,"r")) != NULL, "Can't find %s", path );
- sprintf(path, "%s/snac.plStrain.%d.%06d.restart",context->outputPath,context->rank,context->restartStep);
+ sprintf(path, "%s/snac.plStrain.%d.%06d.restart",context->outputPath,context->rank,context->restartTimestep);
Journal_Firewall( (fp2 = fopen(path,"r")) != NULL, "Can't find %s", path );
/* read in restart file to construct the initial temperature */
for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
Snac_Element* element = Snac_Element_At( context, element_lI );
SnacViscoPlastic_Element* viscoplasticElement = ExtensionManager_Get( context->mesh->elementExtensionMgr, element,
- SnacViscoPlastic_ElementHandle );
+ SnacViscoPlastic_ElementHandle );
const Snac_Material* material = &context->materialProperty[element->material_I];
double plStrain;
- vis_min = context->materialProperty[element->material_I].vis_min;
- if( material->yieldcriterion == druckerprager ) {
+ vis_min = context->materialProperty[element->material_I].vis_min;
+ if( material->yieldcriterion == mohrcoulomb ) {
Tetrahedra_Index tetra_I;
double depls = 0.0f;
for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- double S[3][3];
- Index i,j;
- double depm;
- fscanf( fp1, "%le %le %le %le %le %le %le %le %le", &S[0][0],&S[0][1],&S[0][2],&S[1][0],&S[1][1],&S[1][2],&S[2][0],&S[2][1],&S[2][2]);
- for(i=0;i<3;i++)
- for(j=0;j<3;j++)
- viscoplasticElement->plasticstrainTensor[tetra_I][i][j] = S[i][j];
- /* not the actual viscosity at restartStep, but doesn't affect the calculation afterwards */
- viscoplasticElement->viscosity[tetra_I] = vis_min;
-
- depm = ( S[0][0]+S[1][1]+S[2][2] ) / 3.0f;
- viscoplasticElement->plasticStrain[tetra_I] = sqrt( 0.5f* ((S[0][0]-depm) * (S[0][0]-depm) + (S[1][1]-depm) * (S[1][1]-depm) + (S[2][2]-depm) * (S[2][2]-depm)) + S[0][1]*S[0][1] + S[0][2]*S[0][2] + S[1][2]*S[1][2]);
- } // for
- } // if
- else if( material->yieldcriterion == mohrcoulomb ) {
- Tetrahedra_Index tetra_I;
- double depls = 0.0f;
-
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
double tetraStrain;
Index i,j;
double depm;
fscanf( fp1, "%le", &tetraStrain );
viscoplasticElement->plasticStrain[tetra_I] = tetraStrain;
- /* not the actual viscosity at restartStep, but doesn't affect the calculation afterwards */
+ /* not the actual viscosity at restartTimestep, but doesn't affect the calculation afterwards */
viscoplasticElement->viscosity[tetra_I] = vis_min;
}//for
}//elseif
@@ -165,15 +128,15 @@
for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
Snac_Element* element = Snac_Element_At( context, element_lI );
SnacViscoPlastic_Element* viscoplasticElement = ExtensionManager_Get( context->mesh->elementExtensionMgr, element,
- SnacViscoPlastic_ElementHandle );
- double plStrain = viscoplasticElement->aps;
- vis_min = context->materialProperty[element->material_I].vis_min;
+ SnacViscoPlastic_ElementHandle );
+ double plStrain = viscoplasticElement->aps;
+ vis_min = context->materialProperty[element->material_I].vis_min;
Tetrahedra_Index tetra_I;
for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
viscoplasticElement->plasticStrain[tetra_I] = plStrain;
viscoplasticElement->viscosity[tetra_I] = vis_min;
}//for
}//for
-}//else
+ }//else
}//function
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Output.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Output.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Output.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -42,35 +42,149 @@
#include <math.h>
#include <assert.h>
+
+void _SnacViscoPlastic_WritePlasticStrain( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+
+ if( isTimeToDump( context ) )
+ _SnacViscoPlastic_DumpPlasticStrain( context );
+
+ if( isTimeToCheckpoint( context ) )
+ _SnacViscoPlastic_CheckpointPlasticStrain( context );
+}
+
+
void _SnacViscoPlastic_DumpPlasticStrain( void* _context ) {
Snac_Context* context = (Snac_Context*) _context;
- SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
- context->extensionMgr,
- context,
- SnacViscoPlastic_ContextHandle );
+ SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacViscoPlastic_ContextHandle );
+ Element_LocalIndex element_lI;
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Element* element = Snac_Element_At( context, element_lI );
+ SnacViscoPlastic_Element* elementExt = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ element,
+ SnacViscoPlastic_ElementHandle );
+ float plasticStrain = elementExt->aps;
+ fwrite( &plasticStrain, sizeof(float), 1, contextExt->plStrainOut );
+ }
+ fflush( contextExt->plStrainOut );
+}
- if( context->timeStep ==0 || (context->timeStep-1) % context->dumpEvery == 0 ) {
- Element_LocalIndex element_lI;
- #if DEBUG
- printf( "In %s()\n", __func__ );
- #endif
+void _SnacViscoPlastic_CheckpointPlasticStrain( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacViscoPlastic_ContextHandle );
+ Element_LocalIndex element_lI;
- for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
- Snac_Element* element = Snac_Element_At( context, element_lI );
- SnacViscoPlastic_Element* elementExt = ExtensionManager_Get(
- context->mesh->elementExtensionMgr,
- element,
- SnacViscoPlastic_ElementHandle );
- float plasticStrain = elementExt->aps;
- /* Take average of tetra plastic strain for the element */
- fwrite( &plasticStrain, sizeof(float), 1, contextExt->plStrainOut );
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Element* element = Snac_Element_At( context, element_lI );
+ SnacViscoPlastic_Element* elementExt = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ element,
+ SnacViscoPlastic_ElementHandle );
+ Tetrahedra_Index tetra_I;
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ float plasticStrain = elementExt->plasticStrain[tetra_I];
+ fwrite( &plasticStrain, sizeof(float), 1, contextExt->plStrainCheckpoint );
}
- fflush( contextExt->plStrainOut );
+ fflush( contextExt->plStrainCheckpoint );
}
}
+
+void _SnacViscoPlastic_WriteViscosity( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+
+ if( isTimeToDump( context ) )
+ _SnacViscoPlastic_DumpViscosity( context );
+
+ if( isTimeToCheckpoint( context ) )
+ _SnacViscoPlastic_CheckpointViscosity( context );
+}
+
+
+void _SnacViscoPlastic_DumpViscosity( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacViscoPlastic_ContextHandle );
+ Element_LocalIndex element_lI;
+
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Element* element = Snac_Element_At( context, element_lI );
+ SnacViscoPlastic_Element* elementExt = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ element,
+ SnacViscoPlastic_ElementHandle );
+ Tetrahedra_Index tetra_I;
+ double viscosity = 0.0f;
+ float logviscosity = 0.0f;
+
+ /* Take average of tetra viscosity for the element */
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ )
+ viscosity += elementExt->viscosity[tetra_I]/Tetrahedra_Count;
+ assert(viscosity>0.0);
+ logviscosity = log10(viscosity);
+ fwrite( &logviscosity, sizeof(float), 1, contextExt->viscOut );
+ }
+ fflush( contextExt->viscOut );
+}
+
+
+void _SnacViscoPlastic_CheckpointViscosity( void* _context ) {
+ Snac_Context* context = (Snac_Context*) _context;
+ SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
+ context->extensionMgr,
+ context,
+ SnacViscoPlastic_ContextHandle );
+ Element_LocalIndex element_lI;
+
+#if DEBUG
+ printf( "In %s()\n", __func__ );
+#endif
+
+ for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Element* element = Snac_Element_At( context, element_lI );
+ SnacViscoPlastic_Element* elementExt = ExtensionManager_Get(
+ context->mesh->elementExtensionMgr,
+ element,
+ SnacViscoPlastic_ElementHandle );
+ Tetrahedra_Index tetra_I;
+ double viscosity = 0.0f;
+ float logviscosity = 0.0f;
+
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ viscosity = elementExt->viscosity[tetra_I];
+ assert(viscosity>0.0);
+ logviscosity = log10(viscosity);
+ fwrite( &logviscosity, sizeof(float), 1, contextExt->viscCheckpoint );
+ }
+ }
+ fflush( contextExt->viscCheckpoint );
+}
+
+
+#if 0
void _SnacViscoPlastic_DumpPlasticStrainTensor( void* _context ) {
Snac_Context* context = (Snac_Context*) _context;
SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
@@ -118,38 +232,4 @@
fflush( contextExt->plstrainTensorOut );
}
}
-
-void _SnacViscoPlastic_DumpViscosity( void* _context ) {
- Snac_Context* context = (Snac_Context*) _context;
- SnacViscoPlastic_Context* contextExt = ExtensionManager_Get(
- context->extensionMgr,
- context,
- SnacViscoPlastic_ContextHandle );
-
- if( context->timeStep ==0 || (context->timeStep-1) % context->dumpEvery == 0 ) {
- Element_LocalIndex element_lI;
-
-#if DEBUG
- printf( "In %s()\n", __func__ );
#endif
-
- for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
- Snac_Element* element = Snac_Element_At( context, element_lI );
- SnacViscoPlastic_Element* elementExt = ExtensionManager_Get(
- context->mesh->elementExtensionMgr,
- element,
- SnacViscoPlastic_ElementHandle );
- /* Take average of tetra viscosity for the element */
- Tetrahedra_Index tetra_I;
- double viscosity = 0.0f;
- float logviscosity = 0.0f;
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ )
- viscosity += elementExt->viscosity[tetra_I]/Tetrahedra_Count;
- assert(viscosity>0.0);
- logviscosity = log10(viscosity);
- fwrite( &logviscosity, sizeof(float), 1, contextExt->viscOut );
- }
- fflush( contextExt->viscOut );
- }
-}
-
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Output.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Output.h 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Output.h 2009-04-10 20:43:44 UTC (rev 14656)
@@ -32,8 +32,12 @@
#ifndef _SnacViscoPlastic_Output_
#define _SnacViscoPlastic_Output_
+ void _SnacViscoPlastic_WritePlasticStrain( void* _context );
+ void _SnacViscoPlastic_WriteViscosity( void* _context );
void _SnacViscoPlastic_DumpPlasticStrain( void* _context );
+ void _SnacViscoPlastic_CheckpointPlasticStrain( void* _context );
+ void _SnacViscoPlastic_DumpViscosity( void* _context );
+ void _SnacViscoPlastic_CheckpointViscosity( void* _context );
void _SnacViscoPlastic_DumpPlasticStrainTensor( void* _context );
- void _SnacViscoPlastic_DumpViscosity( void* _context );
#endif
Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Register.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Register.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplastic/Register.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -109,33 +109,23 @@
SnacViscoPlastic_Type );
EntryPoint_Prepend( /* Dump the initial plastic strain */
Context_GetEntryPoint( context, AbstractContext_EP_Execute ),
- "SnacViscoPlastic_Dump",
- _SnacViscoPlastic_DumpPlasticStrain,
+ "SnacViscoPlastic_WritePlasticStrain",
+ _SnacViscoPlastic_WritePlasticStrain,
SnacViscoPlastic_Type );
EntryPoint_Append( /* and dump each loop */
Context_GetEntryPoint( context, Snac_EP_CalcStresses ),
- "SnacViscoPlastic_Dump",
- _SnacViscoPlastic_DumpPlasticStrain,
+ "SnacViscoPlastic_WritePlasticStrain",
+ _SnacViscoPlastic_WritePlasticStrain,
SnacViscoPlastic_Type );
- EntryPoint_Prepend( /* Dump the initial tetra plastic strain */
+ EntryPoint_Prepend( /* Dump the initial viscosity */
Context_GetEntryPoint( context, AbstractContext_EP_Execute ),
- "SnacViscoPlastic_Dump",
- _SnacViscoPlastic_DumpPlasticStrainTensor,
+ "SnacViscoPlastic_WriteViscosity",
+ _SnacViscoPlastic_WriteViscosity,
SnacViscoPlastic_Type );
- EntryPoint_Append( /* and dump tetra plastic strain each loop */
- Context_GetEntryPoint( context, Snac_EP_CalcStresses ),
- "SnacViscoPlastic_Dump",
- _SnacViscoPlastic_DumpPlasticStrainTensor,
- SnacViscoPlastic_Type );
- EntryPoint_Prepend( /* Dump the initial viscosity*/
- Context_GetEntryPoint( context, AbstractContext_EP_Execute ),
- "SnacViscoPlastic_Dump",
- _SnacViscoPlastic_DumpViscosity,
- SnacViscoPlastic_Type );
EntryPoint_Append( /* and dump visocisty each loop */
Context_GetEntryPoint( context, Snac_EP_CalcStresses ),
- "SnacViscoPlastic_Dump",
- _SnacViscoPlastic_DumpViscosity,
+ "SnacViscoPlastic_WriteViscosity",
+ _SnacViscoPlastic_WriteViscosity,
SnacViscoPlastic_Type );
/* Add extensions to the interpolate element entry point, but it will only exist if the remesher is loaded. */
Modified: long/3D/SNAC/trunk/Snac/plugins/winkler/Force.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/winkler/Force.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/winkler/Force.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -149,7 +149,7 @@
(*force)[2] += factor4 * ( press_norm * area * normal[2] );
}
}
- if( context->restartStep == 0 ) {
+ if( context->restartTimestep == 0 ) {
if( context->timeStep == 1 ) {
Fy = (*force)[1];
if(Fy != 0.0)
@@ -563,7 +563,7 @@
}
}
}
- if( context->restartStep == 0 ) {
+ if( context->restartTimestep == 0 ) {
if( context->timeStep == 1 ) {
Fr = (*force)[0]*sin(theta)*cos(phi) + (*force)[1]*sin(theta)*sin(phi) + (*force)[2]*cos(theta);
node->residualFr = Fr;
Modified: long/3D/SNAC/trunk/Snac/plugins/winkler/InitialConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/winkler/InitialConditions.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/winkler/InitialConditions.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -46,7 +46,7 @@
Snac_Context* context = (Snac_Context*)_context;
- if( context->restartStep > 0 && (context->timeStep - context->restartStep == 1) ) {
+ if( context->restartTimestep > 0 && (context->timeStep - context->restartTimestep == 1) ) {
FILE* fp;
char path[PATH_MAX];
Modified: long/3D/SNAC/trunk/Snac/plugins/winkler/Output.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/winkler/Output.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/plugins/winkler/Output.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -48,7 +48,12 @@
Snac_Context* context = (Snac_Context*)_context;
- if( context->timeStep - context->restartStep == 1 ) {
+ /* if( context->timeStep - context->restartTimestep == 1 ) { */
+ /*
+ * It seems fine to record isoForce.* only once at time step = 1.
+ * although there might be a pitfall to be found yet.
+ */
+ if( context->timeStep - context->restartTimestep == 1 ) {
FILE* fp;
char path[PATH_MAX];
Node_LocalIndex node_lI;
Modified: long/3D/SNAC/trunk/Snac/snac2vtk/snac2vtk.c
===================================================================
--- long/3D/SNAC/trunk/Snac/snac2vtk/snac2vtk.c 2009-04-10 18:54:02 UTC (rev 14655)
+++ long/3D/SNAC/trunk/Snac/snac2vtk/snac2vtk.c 2009-04-10 20:43:44 UTC (rev 14656)
@@ -54,7 +54,7 @@
#define PATH_MAX 1024
#endif
#ifndef Tetrahedra_Count
- #define Tetrahedra_Count 10
+ #define Tetrahedra_Count 1
#endif
#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
@@ -124,7 +124,7 @@
int doVisc = 1;
double failureAngle = 30.0;
const unsigned numStressVectorComponent = 6; /* 6 components in stress vector */
-const unsigned numStressComponentsPerElement = 60; /* 6 components times 10 tets per element */
+const unsigned numStressComponentsPerElement = 6; /* 6 averaged components per element */
int main( int argc, char* argv[])
{
@@ -1231,45 +1231,43 @@
double failureAngleRadians=elementStressMeasures->failureAngle*M_PI/180.0;
double normalVector[3],slopeParallelVector[3],tractionVector[3];
double tmp;
+ float stressTensorArray[numStressVectorComponent];
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- float stressTensorArray[numStressVectorComponent];
-
- /*
- * Read all tetrahedral stress tensors for this element gI
- */
- if ( fread( stressTensorArray, sizeof(float), numStressVectorComponent, stressTensorIn )==0 ) {
- if (feof(stressTensorIn)) {
- fprintf(stderr, "Error (reached EOF prematurely) while reading Snac stress tensor output file: tetrahedral element #%d\n" , tetra_I);
- exit(1);
- } else {
- fprintf(stderr, "Error while reading Snac stress tensor output file: tetrahedral element #%d\n" , tetra_I);
- exit(1);
- }
+ /*
+ * Read all tetrahedral stress tensors for this element gI
+ */
+ if ( fread( stressTensorArray, sizeof(float), numStressVectorComponent, stressTensorIn )==0 ) {
+ if (feof(stressTensorIn)) {
+ fprintf(stderr, "Error (reached EOF prematurely) while reading Snac stress tensor output file: tetrahedral element #%d\n" , tetra_I);
+ exit(1);
+ } else {
+ fprintf(stderr, "Error while reading Snac stress tensor output file: tetrahedral element #%d\n" , tetra_I);
+ exit(1);
}
- /*
- * Report error and bail if we pick up NaNs in any of the stress components.
- */
- if(isnan(stressTensorArray[0]) || isnan(stressTensorArray[1])
- || isnan(stressTensorArray[2]) || isnan(stressTensorArray[3])
- || isnan(stressTensorArray[4]) || isnan(stressTensorArray[5])) {
- fprintf(stderr,"NaN in stress tensor file\n");
- abort();
- }
- /*
- * Build average stress tensor for element by summing tetrahedral tensor components
- * - even though it's symmetric, do for all 9 components in case we pick the wrong ones before diagonalization
- */
- elementStressTensor[0][0]+=stressTensorArray[0]/(double)Tetrahedra_Count;
- elementStressTensor[1][1]+=stressTensorArray[1]/(double)Tetrahedra_Count;
- elementStressTensor[2][2]+=stressTensorArray[2]/(double)Tetrahedra_Count;
- elementStressTensor[0][1]+=stressTensorArray[3]/(double)Tetrahedra_Count;
- elementStressTensor[0][2]+=stressTensorArray[4]/(double)Tetrahedra_Count;
- elementStressTensor[1][2]+=stressTensorArray[5]/(double)Tetrahedra_Count;
- elementStressTensor[1][0]+=stressTensorArray[3]/(double)Tetrahedra_Count;
- elementStressTensor[2][0]+=stressTensorArray[4]/(double)Tetrahedra_Count;
- elementStressTensor[2][1]+=stressTensorArray[5]/(double)Tetrahedra_Count;
- }
+ }
+ /*
+ * Report error and bail if we pick up NaNs in any of the stress components.
+ */
+ if(isnan(stressTensorArray[0]) || isnan(stressTensorArray[1])
+ || isnan(stressTensorArray[2]) || isnan(stressTensorArray[3])
+ || isnan(stressTensorArray[4]) || isnan(stressTensorArray[5])) {
+ fprintf(stderr,"NaN in stress tensor file\n");
+ abort();
+ }
+ /*
+ * Build average stress tensor for element by summing tetrahedral tensor components
+ * - even though it's symmetric, do for all 9 components in case we pick the wrong ones before diagonalization
+ */
+ elementStressTensor[0][0]=stressTensorArray[0];
+ elementStressTensor[1][1]=stressTensorArray[1];
+ elementStressTensor[2][2]=stressTensorArray[2];
+ elementStressTensor[0][1]=stressTensorArray[3];
+ elementStressTensor[0][2]=stressTensorArray[4];
+ elementStressTensor[1][2]=stressTensorArray[5];
+ elementStressTensor[1][0]=stressTensorArray[3];
+ elementStressTensor[2][0]=stressTensorArray[4];
+ elementStressTensor[2][1]=stressTensorArray[5];
+
/*
* Diagonalize and find principal stresses from mean stress tensor for element
*/
More information about the CIG-COMMITS
mailing list