[cig-commits] r18875 - long/3D/SNAC/trunk/Snac/libSnac/src
echoi at geodynamics.org
echoi at geodynamics.org
Wed Sep 7 09:28:30 PDT 2011
Author: echoi
Date: 2011-09-07 09:28:30 -0700 (Wed, 07 Sep 2011)
New Revision: 18875
Added:
long/3D/SNAC/trunk/Snac/libSnac/src/Restart.c
long/3D/SNAC/trunk/Snac/libSnac/src/Restart.h
Modified:
long/3D/SNAC/trunk/Snac/libSnac/src/Context.c
long/3D/SNAC/trunk/Snac/libSnac/src/Context.h
long/3D/SNAC/trunk/Snac/libSnac/src/Makefile.def
long/3D/SNAC/trunk/Snac/libSnac/src/Snac.h
long/3D/SNAC/trunk/Snac/libSnac/src/Stress.c
long/3D/SNAC/trunk/Snac/libSnac/src/UpdateElement.c
long/3D/SNAC/trunk/Snac/libSnac/src/UpdateElement.h
Log:
Improved restart functionality:
1. Moved functions reading restarting values from the SnacRestart plugin to the main Snac library: Restart.c/h.
As a result, the SnacRestart plugin is not needed any more and will be removed.
2. Restarting takes a separate route from the normal initialization. The new restarting sequence is implemented in
functions defined in Restart.c,
_Snac_Context_LoopElements_Restart() in Context.c, and
Snac_UpdateElementMomentum_Restart() in UpdateElement.c.
3. New files added to Makefile.def.
4. One trivial change to the sign of hydroPressure in Stress.c.
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Context.c
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Context.c 2011-09-06 02:01:06 UTC (rev 18874)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Context.c 2011-09-07 16:28:30 UTC (rev 18875)
@@ -48,6 +48,7 @@
#include "Force.h"
#include "UpdateNode.h"
#include "Parallel.h"
+#include "Restart.h"
#include "Context.h"
#include <stdio.h>
#include <stdlib.h>
@@ -994,7 +995,7 @@
/* Write out timeStep info */
_Snac_Context_WriteLoopInfo( self );
- if( self->restartTimestep == 0 ) {
+ if( self->restartTimestep == 0 ) { /* If not restarting (i.e., starting normally) */
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;
@@ -1017,11 +1018,17 @@
}
element->stress = 0.5f * sqrt( 0.5f * fabs( sVolAvg/Tetrahedra_Count ) );
}
+ /* 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 );
}
+ else if( (self->restartTimestep > 0) && (self->timeStep==self->restartTimestep) ) { /* if restarting */
+ _Snac_Restart_ResetMinLengthScale( self );
+ _Snac_Restart_InitialCoords( self );
+ _Snac_Restart_InitialVelocities( self );
+ _Snac_Context_LoopElements_Restart( self );
+ _Snac_Restart_InitialStress( 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 );
-
_Snac_Context_WriteOutput( self );
KeyCall( self, self->syncK, EntryPoint_Class_VoidPtr_CallCast* )( KeyHandle(self,self->syncK), self );
@@ -1176,7 +1183,28 @@
}
}
+void _Snac_Context_LoopElements_Restart( void* context ) {
+ Snac_Context* self = (Snac_Context*)context;
+ Element_LocalIndex element_lI;
+ if( self->rank == 0 ) Journal_DPrintf( self->debug, "In: %s\n", __func__ );
+ if( self->rank == 0 ) Journal_Printf(
+ self->verbose,
+ "For each element, calc volume, surface vel, and then min length scale\n" );
+
+ /* Update all the elements, and in the process work out this processor's minLengthScale */
+ element_lI = 0;
+ Snac_UpdateElementMomentum_Restart( self, element_lI, &self->minLengthScale );
+ for( element_lI = 1; element_lI < self->mesh->elementLocalCount; element_lI++ ) {
+ double elementMinLengthScale;
+ Snac_UpdateElementMomentum_Restart( self, element_lI, &elementMinLengthScale );
+ if( elementMinLengthScale < self->minLengthScale ) {
+ self->minLengthScale = elementMinLengthScale;
+ }
+ }
+}
+
+
void _Snac_Context_Sync( void* context ) {
Snac_Context* self = (Snac_Context*)context;
Node_DomainIndex node_dI;
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Context.h 2011-09-06 02:01:06 UTC (rev 18874)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Context.h 2011-09-07 16:28:30 UTC (rev 18875)
@@ -224,6 +224,7 @@
/* Solve volumes, normals, etc for each element */
void _Snac_Context_LoopElements( void* context );
+ void _Snac_Context_LoopElements_Restart( void* context );
/* Sync for the Snac implementation */
void _Snac_Context_Sync( void* context );
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Makefile.def
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Makefile.def 2011-09-06 02:01:06 UTC (rev 18874)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Makefile.def 2011-09-07 16:28:30 UTC (rev 18875)
@@ -46,6 +46,7 @@
Force.c \
UpdateNode.c \
Parallel.c \
+ Restart.c \
Context.c \
SnacSync.c \
Init.c \
@@ -63,6 +64,7 @@
Element.h \
EntryPoint.h \
UpdateElement.h \
+ Restart.h \
StrainRate.h \
Stress.h \
Force.h \
Added: long/3D/SNAC/trunk/Snac/libSnac/src/Restart.c
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Restart.c (rev 0)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Restart.c 2011-09-07 16:28:30 UTC (rev 18875)
@@ -0,0 +1,250 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+** $Id: InitialConditions.c 2498 2005-01-07 03:57:00Z PatrickSunter $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StGermain/FD/FD.h>
+#include "Snac/Snac.h"
+#include "Restart.h"
+
+#include <string.h>
+#ifndef PATH_MAX
+ #define PATH_MAX 1024
+#endif
+
+#define DEBUG
+
+void _Snac_Restart_ResetMinLengthScale( void* _context ) {
+
+ Snac_Context* context = (Snac_Context*)_context;
+ FILE* fp;
+ char fname[PATH_MAX];
+ sprintf( fname, "%s/snac.minLengthScale.restart", context->outputPath );
+#ifdef DEBUG
+ fprintf(stderr,"Restarter: reading min length scale from %s\n",fname);
+#endif
+ Journal_Firewall( (fp = fopen( fname, "r" )) != NULL, "Failed to open %s\n", fname );
+ fscanf(fp,"%le",&(context->initMinLengthScale));
+ fclose( fp );
+
+}
+
+void _Snac_Restart_InitialCoords( void* _context ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ FILE* fp;
+ Node_LocalIndex node_lI;
+ char path[PATH_MAX];
+
+ Journal_Printf( context->snacInfo, "In: %s\n", __func__ );
+
+ 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->restartTimestep);
+#endif
+ 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++ ) {
+ Coord* coord = Snac_NodeCoord_P( context, node_lI );
+ double x,y,z;
+ fscanf( fp, "%le %le %le", &x,&y,&z);
+
+ (*coord)[0] = x;
+ (*coord)[1] = y;
+ (*coord)[2] = z;
+#ifdef DEBUG
+/* fprintf(stderr,"Node %d: %g, %g, %g\n", node_lI, x,y,z); */
+#endif
+ }
+ fclose( fp );
+}
+
+void _Snac_Restart_InitialVelocities( void* _context ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ FILE* fp;
+ Node_LocalIndex node_lI;
+ char path[PATH_MAX];
+
+ Journal_Printf( context->snacInfo, "In: %s\n", __func__ );
+
+ 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->restartTimestep);
+#endif
+ 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++ ) {
+ Snac_Node* node = Snac_Node_At( context, node_lI );
+ double vx,vy,vz;
+ fscanf( fp, "%le %le %le", &vx,&vy,&vz);
+
+ node->velocity[0] = vx;
+ node->velocity[1] = vy;
+ node->velocity[2] = vz;
+ }
+ fclose( fp );
+}
+
+#if 0
+/* void _Snac_Restart_InitialStress( void* _context ) { */
+/* Snac_Context* context = (Snac_Context*)_context; */
+/* FILE* fp; */
+/* Element_LocalIndex element_lI; */
+/* 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->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->restartTimestep); */
+/* #endif */
+
+/* /\* read in restart file to assign initial stress field. *\/ */
+/* for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) { */
+/* Snac_Element* element = Snac_Element_At( context, element_lI ); */
+/* Tetrahedra_Index tetra_I; */
+/* Stress sVolAvg=0.0f; */
+/* Stress sOtherAvg=0.0f; */
+
+/* for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) { */
+/* double S[3][3]; */
+/* Index i,j; */
+/* fscanf( fp, "%le %le %le %le %le %le", &S[0][0],&S[1][1],&S[2][2],&S[0][1],&S[0][2],&S[1][2]); */
+/* S[1][0] = S[0][1]; */
+/* S[2][0] = S[0][2]; */
+/* S[2][1] = S[1][2]; */
+/* for(i=0;i<3;i++) */
+/* for(j=0;j<3;j++) */
+/* element->tetra[tetra_I].stress[i][j] = S[i][j]; */
+/* } */
+/* element->hydroPressure = 0.0; */
+/* for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) { */
+/* 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] + */
+/* element->tetra[tetra_I].stress[0][0] * element->tetra[tetra_I].stress[1][1]; */
+/* sOtherAvg += */
+/* element->tetra[tetra_I].stress[0][1] * element->tetra[tetra_I].stress[0][1] + */
+/* element->tetra[tetra_I].stress[1][2] * element->tetra[tetra_I].stress[1][2] + */
+/* element->tetra[tetra_I].stress[0][2] * element->tetra[tetra_I].stress[0][2]; */
+/* } */
+/* sVolAvg /= Tetrahedra_Count; */
+/* sOtherAvg /= Tetrahedra_Count; */
+/* element->stress = 0.5 * sqrt( 0.5 * fabs( -1.0f * sVolAvg + sOtherAvg ) ); */
+/* element->hydroPressure *= -1.0; /\* for consistent sign convention. *\/ */
+/* } */
+/* fclose( fp ); */
+/* } */
+#endif
+
+void _Snac_Restart_InitialStress( void* _context ) {
+ Snac_Context* context = (Snac_Context*)_context;
+ FILE* fp;
+ Element_LocalIndex element_lI;
+ 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->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->restartTimestep);
+#endif
+
+ /* read in restart file to assign initial stress field. */
+ for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
+ Snac_Element* element = Snac_Element_At( context, element_lI );
+ Tetrahedra_Index tetra_I;
+ Stress traceStress[Tetrahedra_Count];
+ Stress partialStress;
+ Stress sVolAvg=0.0f;
+ Stress sOtherAvg=0.0f;
+ Stress pressure=0.0f;
+ double elemVolume = ( Tetrahedra_Count > 5 )?(2.0*element->volume):element->volume;
+
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ double S[3][3];
+ Index i,j;
+ fscanf( fp, "%le %le %le %le %le %le", &S[0][0],&S[1][1],&S[2][2],&S[0][1],&S[0][2],&S[1][2]);
+ S[1][0] = S[0][1];
+ S[2][0] = S[0][2];
+ S[2][1] = S[1][2];
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++)
+ element->tetra[tetra_I].stress[i][j] = S[i][j];
+ }
+
+
+ partialStress = 0.0f;
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+
+ /* Initialise the trace stress tensor for this tetrahedra to 0 */
+ memset( &traceStress[tetra_I], 0, sizeof(Stress) );
+ traceStress[tetra_I]= (
+ element->tetra[tetra_I].stress[0][0] +
+ element->tetra[tetra_I].stress[1][1] +
+ element->tetra[tetra_I].stress[2][2] ) / 3.0f;
+ partialStress += traceStress[tetra_I] * element->tetra[tetra_I].volume;
+ }
+
+ sVolAvg = 0.0f;
+ sOtherAvg = 0.0f;
+ pressure = 0.0f;
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ const double averageTraceStress = partialStress / elemVolume;
+ element->tetra[tetra_I].stress[0][0] += -1.0f * traceStress[tetra_I] + averageTraceStress;
+ element->tetra[tetra_I].stress[1][1] += -1.0f * traceStress[tetra_I] + averageTraceStress;
+ element->tetra[tetra_I].stress[2][2] += -1.0f * traceStress[tetra_I] + averageTraceStress;
+
+ sVolAvg +=
+ (element->tetra[tetra_I].stress[1][1]-averageTraceStress) * (element->tetra[tetra_I].stress[2][2]-averageTraceStress) +
+ (element->tetra[tetra_I].stress[2][2]-averageTraceStress) * (element->tetra[tetra_I].stress[0][0]-averageTraceStress) +
+ (element->tetra[tetra_I].stress[0][0]-averageTraceStress) * (element->tetra[tetra_I].stress[1][1]-averageTraceStress);
+ sOtherAvg +=
+ element->tetra[tetra_I].stress[0][1] * element->tetra[tetra_I].stress[0][1] +
+ element->tetra[tetra_I].stress[1][2] * element->tetra[tetra_I].stress[1][2] +
+ element->tetra[tetra_I].stress[0][2] * element->tetra[tetra_I].stress[0][2];
+
+ pressure += (element->tetra[tetra_I].stress[0][0]+element->tetra[tetra_I].stress[1][1]+element->tetra[tetra_I].stress[2][2])/3.0f;
+ }
+ sVolAvg /= Tetrahedra_Count;
+ sOtherAvg /= Tetrahedra_Count;
+ pressure /= Tetrahedra_Count;
+
+ /* Calculate the element stress from the tetrahedra stress tensors */
+ element->stress = 0.5f * sqrt( 0.5f * fabs( -1.0f * sVolAvg + sOtherAvg ) );
+ element->hydroPressure = -1.0*pressure;
+ }
+ fclose( fp );
+}
Added: long/3D/SNAC/trunk/Snac/libSnac/src/Restart.h
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Restart.h (rev 0)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Restart.h 2011-09-07 16:28:30 UTC (rev 18875)
@@ -0,0 +1,50 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003,
+** Steve Quenette, 110 Victoria Street, Melbourne, Victoria, 3053, Australia.
+** Californian Institute of Technology, 1200 East California Boulevard, Pasadena, California, 91125, USA.
+** University of Texas, 1 University Station, Austin, Texas, 78712, USA.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Stevan M. Quenette, Visitor in Geophysics, Caltech.
+** Luc Lavier, Research Scientist, The University of Texas. (luc at utig.ug.utexas.edu)
+** Luc Lavier, Research Scientist, Caltech.
+**
+** This program is free software; you can redistribute it and/or modify it
+** under the terms of the GNU General Public License as published by the
+** Free Software Foundation; either version 2, or (at your option) any
+** later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+**
+*/
+/** \file
+** Role:
+** Read and assign values to variables that are necessary for restarting.
+**
+** Assumptions:
+** In general, when restarting, context->restartTimestep > 0 and context->timeStep == restartTimestep.
+**
+** Comments:
+**
+** $Id: Restart.h 1431 2004-05-18 07:19:21Z SteveQuenette $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __Snac_Restart_h__
+#define __Snac_Restart_h__
+
+ void _Snac_Restart_ResetMinLengthScale( void* _context );
+ void _Snac_Restart_InitialCoords( void* _context );
+ void _Snac_Restart_InitialVelocities( void* _context );
+ void _Snac_Restart_InitialStress( void* _context );
+
+#endif /* __Snac_Restart_h__ */
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Snac.h
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Snac.h 2011-09-06 02:01:06 UTC (rev 18874)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Snac.h 2011-09-07 16:28:30 UTC (rev 18875)
@@ -54,6 +54,7 @@
#include "Element.h"
#include "EntryPoint.h"
#include "UpdateElement.h"
+ #include "Restart.h"
#include "StrainRate.h"
#include "Stress.h"
#include "Force.h"
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Stress.c
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Stress.c 2011-09-06 02:01:06 UTC (rev 18874)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Stress.c 2011-09-07 16:28:30 UTC (rev 18875)
@@ -145,7 +145,7 @@
element->tetra[tetra_I].stress[1][2] * element->tetra[tetra_I].stress[1][2] +
element->tetra[tetra_I].stress[0][2] * element->tetra[tetra_I].stress[0][2];
- pressure += -1.0f * (element->tetra[tetra_I].stress[0][0]+element->tetra[tetra_I].stress[1][1]+element->tetra[tetra_I].stress[2][2])/3.0f;
+ pressure += (element->tetra[tetra_I].stress[0][0]+element->tetra[tetra_I].stress[1][1]+element->tetra[tetra_I].stress[2][2])/3.0f;
}
sVolAvg /= Tetrahedra_Count;
sOtherAvg /= Tetrahedra_Count;
@@ -153,7 +153,7 @@
/* Calculate the element stress from the tetrahedra stress tensors */
element->stress = 0.5f * sqrt( 0.5f * fabs( -1.0f * sVolAvg + sOtherAvg ) );
- element->hydroPressure = pressure;
+ element->hydroPressure = -1.0*pressure;
/* To catch nan or inf in the stress values even in the optimised mode. */
/* In the usual DEBUG mode, the above Journal_DFirewalls are sufficient. */
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/UpdateElement.c
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/UpdateElement.c 2011-09-06 02:01:06 UTC (rev 18874)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/UpdateElement.c 2011-09-07 16:28:30 UTC (rev 18875)
@@ -75,7 +75,7 @@
/* Calculate the volume, surface { area, normal and velocity }, and smallest length scale for each tetrahedra. */
for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
Tetrahedra_Surface_Index surface_I;
- double lengthScale;
+ double lengthScale;
/*ccccc*/
memset( element->tetra[tetra_I].strain, 0, sizeof(element->tetra[tetra_I].strain) );
@@ -83,10 +83,10 @@
/* Calculate this tetrahedra's volume */
element->tetra[tetra_I].old_volume = element->tetra[tetra_I].volume;
element->tetra[tetra_I].volume = Tetrahedra_Volume(
- Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][0] ),
- Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][1] ),
- Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][2] ),
- Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][3] ) );
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][0] ),
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][1] ),
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][2] ),
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][3] ) );
if(self->timeStep == 0)
element->tetra[tetra_I].old_volume = element->tetra[tetra_I].volume;
@@ -239,7 +239,7 @@
element->tetra[tetra_I].strain[1][2] += ( rotation[tetra_I][1][2] *
(element->tetra[tetra_I].strain[2][2] - element->tetra[tetra_I].strain[1][1] ) ) * self->dt +
( (-1.0f*rotation[tetra_I][0][1]) * element->tetra[tetra_I].strain[0][2] -
- rotation[tetra_I][0][2] * element->tetra[tetra_I].strain[0][1] ) * self->dt;
+ rotation[tetra_I][0][2] * element->tetra[tetra_I].strain[0][1] ) * self->dt;
/* Rotate the stress of this tetrahedra. */
element->tetra[tetra_I].stress[0][0] += (
@@ -267,3 +267,94 @@
if( Tetrahedra_Count > 5 )
element->volume *= 0.5;
}
+
+void Snac_UpdateElementMomentum_Restart( void* context, Element_LocalIndex element_lI, double* elementMinLengthScale ) {
+ Snac_Context* self = (Snac_Context*)context;
+ Tetrahedra_Index tetra_I;
+ Snac_Element* element = Snac_Element_At( self, element_lI );
+
+ Snac_Material* material = &self->materialProperty[element->material_I];
+
+
+ /* Find this processor's maximum Vp and store it as speedOfSound. */
+ if( self->dtType == Snac_DtType_Wave ) {
+ double Vp = sqrt((material->lambda+2.0f*material->mu)/material->phsDensity);
+ if( Vp > self->speedOfSound )
+ self->speedOfSound = Vp;
+ }
+
+ /* Initialise output data */
+ element->volume=0.0;
+ *elementMinLengthScale = Tetrahedra_Max_Propagation_Length;
+
+ /* Calculate the volume, surface { area, normal and velocity }, and smallest length scale for each tetrahedra. */
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ Tetrahedra_Surface_Index surface_I;
+ double lengthScale;
+
+ /*ccccc*/
+ memset( element->tetra[tetra_I].strain, 0, sizeof(element->tetra[tetra_I].strain) );
+
+ /* Calculate this tetrahedra's volume */
+ element->tetra[tetra_I].old_volume = element->tetra[tetra_I].volume;
+ element->tetra[tetra_I].volume = Tetrahedra_Volume(
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][0] ),
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][1] ),
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][2] ),
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][3] ) );
+ if(self->timeStep == 0)
+ element->tetra[tetra_I].old_volume = element->tetra[tetra_I].volume;
+
+ /* Calculate this tetrahedra's area, normal and velocity for each surface. */
+ for( surface_I = 0; surface_I < Tetrahedra_Surface_Count; surface_I++ ) {
+ element->tetra[tetra_I].surface[surface_I].area = Tetrahedra_SurfaceArea(
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][FaceToNode[surface_I][0]] ),
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][FaceToNode[surface_I][1]] ),
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][FaceToNode[surface_I][2]] ) );
+
+ Tetrahedra_SurfaceNormal(
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][FaceToNode[surface_I][0]] ),
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][FaceToNode[surface_I][1]] ),
+ Snac_Element_NodeCoord( self, element_lI, TetraToNode[tetra_I][FaceToNode[surface_I][2]] ),
+ &element->tetra[tetra_I].surface[surface_I].normal );
+ }
+
+
+ /* keep track of this processor's smallest length scale */
+ lengthScale = element->tetra[tetra_I].surface[0].area;
+ if( lengthScale ) {
+ lengthScale = fabs( element->tetra[tetra_I].volume / lengthScale );
+ if( *elementMinLengthScale > lengthScale ) {
+ *elementMinLengthScale = lengthScale;
+ }
+ }
+
+ lengthScale = element->tetra[tetra_I].surface[1].area;
+ if( lengthScale ) {
+ lengthScale = fabs( element->tetra[tetra_I].volume / lengthScale );
+ if( *elementMinLengthScale > lengthScale ) {
+ *elementMinLengthScale = lengthScale;
+ }
+ }
+
+ lengthScale = element->tetra[tetra_I].surface[2].area;
+ if( lengthScale ) {
+ lengthScale = fabs( element->tetra[tetra_I].volume / lengthScale );
+ if( *elementMinLengthScale > lengthScale ) {
+ *elementMinLengthScale = lengthScale;
+ }
+ }
+
+ lengthScale = element->tetra[tetra_I].surface[3].area;
+ if( lengthScale ) {
+ lengthScale = fabs( element->tetra[tetra_I].volume / lengthScale );
+ if( *elementMinLengthScale > lengthScale ) {
+ *elementMinLengthScale = lengthScale;
+ }
+ }
+ /* Ammend the element volume with this tetrahedra's contribution */
+ element->volume += element->tetra[tetra_I].volume;
+ }
+ if( Tetrahedra_Count > 5 )
+ element->volume *= 0.5;
+}
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/UpdateElement.h
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/UpdateElement.h 2011-09-06 02:01:06 UTC (rev 18874)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/UpdateElement.h 2011-09-07 16:28:30 UTC (rev 18875)
@@ -44,5 +44,6 @@
#define __Snac_UpdateElement_h__
void Snac_UpdateElementMomentum( void* context, Element_LocalIndex element_lI, double* elementMinLengthScale );
+ void Snac_UpdateElementMomentum_Restart( void* context, Element_LocalIndex element_lI, double* elementMinLengthScale );
#endif /* __Snac_UpdateElement_h__ */
More information about the CIG-COMMITS
mailing list