[cig-commits] r16146 - long/3D/SNAC/trunk/Snac/libSnac/src
echoi at geodynamics.org
echoi at geodynamics.org
Tue Jan 19 06:15:51 PST 2010
Author: echoi
Date: 2010-01-19 06:15:51 -0800 (Tue, 19 Jan 2010)
New Revision: 16146
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/Force.c
long/3D/SNAC/trunk/Snac/libSnac/src/UpdateElement.c
Log:
Added an option for "dtType" for mass-scaling and dt calculation.
- in addition to "dynamic" and "constant",
"wave" is a new value for the input key, "dtType".
- The "wave" type make Snac uses real mass without scaling and dt bounded by the global maximum of P-wave velocity.
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Context.c
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Context.c 2010-01-19 14:05:43 UTC (rev 16145)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Context.c 2010-01-19 14:15:51 UTC (rev 16146)
@@ -83,6 +83,7 @@
/* Snac_Context dt type names */
const Name Snac_DtType_Constant = "constant";
const Name Snac_DtType_Dynamic = "dynamic";
+const Name Snac_DtType_Wave = "wave";
Snac_Context* Snac_Context_New(
@@ -321,6 +322,17 @@
Journal_Printf( self->info, "\"dtType\" set by Dictionary to \"%s\"\n", self->dtType );
}
+ else if( !strcmp( tmpStr, Snac_DtType_Wave ) ) {
+ /* Set dtType to "dynamic" using known pointer/key... it means checks are pointer comparisons only, not strcmp */
+ self->dtType = Snac_DtType_Wave;
+
+ /* When Courant Dt type is specified, Dt will be calculated later,
+ upper-bounded by the Courant-Fredrichs-Levy condition. Give 0 as starting value... legacy behaviour. */
+ self->dt = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( self->dictionary, "timeStep", Dictionary_Entry_Value_FromDouble( 0.0 ) ) );
+
+ Journal_Printf( self->info, "\"dtType\" set by Dictionary to \"%s\"\n", self->dtType );
+ }
else {
/* Set dtType to "dynamic" using known pointer/key... it means checks are pointer comparisons only, not strcmp */
/* Warn the user that this assumption has occured */
@@ -1172,6 +1184,11 @@
self->speedOfSound = self->minLengthScale * 0.5f / self->dt;
}
+ else if( self->dtType == Snac_DtType_Wave ) {
+ MPI_Allreduce( &self->speedOfSound, &tmp, 1, MPI_DOUBLE, MPI_MIN, self->communicator );
+ self->speedOfSound = tmp;
+ self->dt = 0.5f * self->minLengthScale / self->speedOfSound;
+ }
if( self->rank == 0 ) {
Journal_Printf( self->debug, "new self->minlengthscale: %g\n", self->minLengthScale );
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Context.h 2010-01-19 14:05:43 UTC (rev 16145)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Context.h 2010-01-19 14:15:51 UTC (rev 16146)
@@ -62,6 +62,7 @@
/* Snac_Context dt type names */
extern const Name Snac_DtType_Constant;
extern const Name Snac_DtType_Dynamic;
+ extern const Name Snac_DtType_Wave;
/* Snac_Context info */
#define __Snac_Context \
@@ -79,7 +80,7 @@
double topoGradMax; \
double topoGradCriterion; \
\
- Snac_Material* materialProperty; \
+ Snac_Material* materialProperty; \
double dt; \
Type dtType; \
double strain_inert; \
@@ -91,9 +92,9 @@
double gravity; \
double demf; \
\
- CompositeVC* nodeICs; \
- CompositeVC* elementICs; \
- CompositeVC* velocityBCs; \
+ CompositeVC* nodeICs; \
+ CompositeVC* elementICs; \
+ CompositeVC* velocityBCs; \
\
/* TODO... we want journal or the like to look after this in the end */ \
FILE* simInfo; \
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Force.c
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Force.c 2010-01-19 14:05:43 UTC (rev 16145)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Force.c 2010-01-19 14:15:51 UTC (rev 16146)
@@ -65,11 +65,7 @@
const double factor4 = 1.0f / 4.0f;
double dirnorm;
double dir2centr[3];
- Index i;
const double gravity = self->gravity;
- int Spherical = 0;
- Dictionary_Entry_Value* pluginsList;
- Dictionary_Entry_Value* plugin;
if(self->spherical) {
Coord* coord = Snac_NodeCoord_P( self, node_lI );
@@ -87,70 +83,43 @@
memset( *force, 0, sizeof(Force) );
memset( *balance, 0, sizeof(Force) );
if( self->forceCalcType == Snac_Force_Complete ) {
- /* Work out the inertial mass contributions */
- if( self->dtType != Snac_DtType_Constant ) {
- /* for each element of this node */
- for( nodeElement_I = 0; nodeElement_I < nodeElementCount; nodeElement_I++ ) {
- Element_DomainIndex element_lI = self->mesh->nodeElementTbl[node_lI][nodeElement_I];
- /* Nodes on the meshes boundary do not have all 8 neighbouring elements... however we are now
- using a shadow depth of 1, so they are in the domain tables */
- if( element_lI < self->mesh->elementDomainCount ) {
- Snac_Element* element = Snac_Element_At( self, element_lI );
- Material_Index material_I = element->material_I;
- const Snac_Material* material = &self->materialProperty[material_I];
- const Density inertialDensity = (material->lambda + 2.0f * material->mu )/
- (speedOfSound * speedOfSound);
- Index index;
-
- Journal_DFirewall(
- inertialDensity > 0.0f,
- self->snacError,
- "forceCalc: Quick, dtType != Constant, element_lI: %u, inertialDensity is less than 0", element_lI );
-
- /* for each tetra that refers to this node */
- for( index = 0; index < Node_Element_Tetrahedra_Count;index++ ) {
- Tetrahedra_Index tetra_I = NodeToTetra[nodeElement_I][index];
-
- *inertialMass += factor4 * inertialDensity * element->tetra[tetra_I].volume;
- Journal_DFirewall(
- !isnan( *inertialMass ) && !isinf( *inertialMass ),
- self->snacError,
- "forceCalc: Quick, dtType: !Constant, element_lI: %u, inertialDensity is either nan or inf", element_lI );
- }
- }
- }
- }
-
for( nodeElement_I = 0; nodeElement_I < nodeElementCount; nodeElement_I++ ) { /* for each element of this node */
Element_DomainIndex element_lI = self->mesh->nodeElementTbl[node_lI][nodeElement_I];
/* Nodes on the meshes boundary do not have all 8 neighbouring elements... but they are stored in
in the domain */
if( element_lI < self->mesh->elementDomainCount ) {
- Snac_Element* element = Snac_Element_At( self, element_lI );
+ Snac_Element* element = Snac_Element_At( self, element_lI );
Material_Index material_I = element->material_I;
- const Snac_Material* material = &self->materialProperty[material_I];
- Index index;
+ const Snac_Material* material = &self->materialProperty[material_I];
+ const Density inertialDensity = (material->lambda + 2.0f * material->mu )/(speedOfSound * speedOfSound);
+ Index index;
/* for each tetra that refers to this node */
for( index = 0; index < Node_Element_Tetrahedra_Count;index++ ) {
- const Tetrahedra_Index tetra_I = NodeToTetra[nodeElement_I][index];
+ const Tetrahedra_Index tetra_I = NodeToTetra[nodeElement_I][index];
const Tetrahedra_Surface_Index surface_I = NodeToSurface[nodeElement_I][index];
/* Element info shortcuts */
- const Snac_Element_Tetrahedra* tetra = &element->tetra[tetra_I];
- const StressTensor* stress = (const StressTensor*) &tetra->stress;
+ const Snac_Element_Tetrahedra* tetra = &element->tetra[tetra_I];
+ const StressTensor* stress = (const StressTensor*) &tetra->stress;
const Snac_Element_Tetrahedra_Surface* surface = &tetra->surface[surface_I];
- const Normal* normal = &surface->normal;
- const Density effDensity = tetra->density;
+ const Normal* normal = &surface->normal;
+ const Density effDensity = tetra->density;
- /*ccccc*/
- if( self->dtType == Snac_DtType_Constant ) {
- double alpha1 = (material->lambda + 2.0f / 3.0f *
- material->mu ) + 4.0f / 3.0f *
- material->mu;
- int dim;
+ /* Work out the mass contributions */
+ if( self->dtType == Snac_DtType_Dynamic ) {
+
+ *inertialMass += factor4 * inertialDensity * element->tetra[tetra_I].volume;
+ Journal_DFirewall(
+ !isnan( *inertialMass ) && !isinf( *inertialMass ),
+ self->snacError,
+ "forceCalc: Complete, dtType: Dynamic, element_lI: %u, inertialDensity is either nan or inf", element_lI );
+ }
+ else if( self->dtType == Snac_DtType_Constant ) {
+ double alpha1 = (material->lambda + 2.0f*material->mu);
+ int dim;
double temp;
double area_sum = 0.0;
@@ -164,8 +133,21 @@
Journal_DFirewall(
!isnan( *inertialMass ) && !isinf( *inertialMass ),
self->snacError,
- "forceCalc: Quick, dtType: Constant, element_lI: %u, inertialDensity is either nan or inf", element_lI );
+ "forceCalc: Complete, dtType: Constant, element_lI: %u, inertialDensity is either nan or inf", element_lI );
}
+ else if( self->dtType == Snac_DtType_Wave ) {
+ *inertialMass += factor4 * effDensity * element->tetra[tetra_I].volume;
+ Journal_DFirewall(
+ !isnan( *inertialMass ) && !isinf( *inertialMass ),
+ self->snacError,
+ "forceCalc: Complete, dtType: Courant, element_lI: %u, inertialDensity is either nan or inf", element_lI );
+ }
+ else {
+ Journal_DFirewall(
+ 0,
+ self->snacError,
+ "forceCalc: Complete, dtType: Unknown. Aborting.\n");
+ }
/*ccccc*/
/* Incorporate this surface's stresses into the force */
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/UpdateElement.c
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/UpdateElement.c 2010-01-19 14:05:43 UTC (rev 16145)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/UpdateElement.c 2010-01-19 14:15:51 UTC (rev 16146)
@@ -57,7 +57,15 @@
Snac_Element* element = Snac_Element_At( self, element_lI );
Rotation rotation[Tetrahedra_Count];
+ 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;
More information about the CIG-COMMITS
mailing list