[cig-commits] r14548 - in long/3D/SNAC/trunk/Snac/plugins: hillSlope plastic
cstark at geodynamics.org
cstark at geodynamics.org
Mon Mar 30 12:26:03 PDT 2009
Author: cstark
Date: 2009-03-30 12:26:02 -0700 (Mon, 30 Mar 2009)
New Revision: 14548
Modified:
long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h
long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.c
long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.c
long/3D/SNAC/trunk/Snac/plugins/plastic/Constitutive.c
Log:
Extended hillslope plugin to handle (1) hetergeneous cohesion and (2) adjustable subdomain position. Also (3) modified plugin registration of CreateWeakPoints.c so that it has an entry point that ensures weak point seeding before stress/strain computation in that time step.
Details:
1) Added input.xml variables "subDomainMinCohesion" and "subDomainMaxCohesion". If these are non-negative the subdomain is filled with random variates spanning this range. The sampling pdf is triangular generated by summing two uniform variates produced by drand48 (ick).
To turn this off, and simultaneously turn off weak point seeding, make "fractionWeakPoints" negative.
To turn this on, but turn off weak point seeding, make "fractionWeakPoints" zero.
To this off but keep weak point seeding on, make "fractionWeakPoints" non-negative (positive inc zero) AND make "subDomainMinCohesion" negative.
2) Added input.xml variable "xSubDomainLeft", which is itself a fraction giving the relative position of the left edge of the subdomain.
Basic error checking is done to ensure that the subdomain can't start too far left (<0) or spill over beyond the domain right edge.
3) In Register.c shifted entry point of CreateWeakPoints.c to Snac_EP_Constitutive. Has the effect of calling CreateWeakPoints.c during the update phase ahead of stress/strain computations - without changing the plugin order in the input.xml.
Practical consequence is that a dump of the first time step at which CreateWeakPoints is run generates a plastic strain output file that reflects its effect.
Modified: long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c 2009-03-30 15:39:22 UTC (rev 14547)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/ConstructExtensions.c 2009-03-30 19:26:02 UTC (rev 14548)
@@ -107,31 +107,42 @@
contextExt->xSubDomainFraction = Dictionary_Entry_Value_AsDouble(
Dictionary_GetDefault( context->dictionary, "xSubDomainFraction",
- Dictionary_Entry_Value_FromDouble( 0.5f ) ) );
+ Dictionary_Entry_Value_FromDouble( 0.7f ) ) );
contextExt->ySubDomainFraction = Dictionary_Entry_Value_AsDouble(
Dictionary_GetDefault( context->dictionary, "ySubDomainFraction",
- Dictionary_Entry_Value_FromDouble( 0.4f ) ) );
+ Dictionary_Entry_Value_FromDouble( 0.8f ) ) );
contextExt->zSubDomainFraction = Dictionary_Entry_Value_AsDouble(
Dictionary_GetDefault( context->dictionary, "zSubDomainFraction",
- Dictionary_Entry_Value_FromDouble( 0.75f ) ) );
+ Dictionary_Entry_Value_FromDouble( 1.0f ) ) );
+ contextExt->xSubDomainLeft = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "xSubDomainLeft",
+ Dictionary_Entry_Value_FromDouble( -0.15f ) ) );
+
contextExt->xTriggerPointFraction = Dictionary_Entry_Value_AsDouble(
Dictionary_GetDefault( context->dictionary, "xTriggerPointFraction",
Dictionary_Entry_Value_FromDouble( 0.5f ) ) );
contextExt->yTriggerPointFraction = Dictionary_Entry_Value_AsDouble(
Dictionary_GetDefault( context->dictionary, "yTriggerPointFraction",
- Dictionary_Entry_Value_FromDouble( 0.4f ) ) );
+ Dictionary_Entry_Value_FromDouble( 0.3f ) ) );
contextExt->zTriggerPointFraction = Dictionary_Entry_Value_AsDouble(
Dictionary_GetDefault( context->dictionary, "zTriggerPointFraction",
Dictionary_Entry_Value_FromDouble( 0.5f ) ) );
contextExt->weakPointCohesion = Dictionary_Entry_Value_AsDouble(
Dictionary_GetDefault( context->dictionary, "weakPointCohesion",
- Dictionary_Entry_Value_FromDouble( 4.0e+05f ) ) );
+ Dictionary_Entry_Value_FromDouble( 2.0e+04f ) ) );
contextExt->triggerPointCohesion = Dictionary_Entry_Value_AsDouble(
Dictionary_GetDefault( context->dictionary, "triggerPointCohesion",
- Dictionary_Entry_Value_FromDouble( 1.0e+04f ) ) );
+ Dictionary_Entry_Value_FromDouble( 2.0e+04f ) ) );
+ contextExt->subDomainMinCohesion = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "subDomainMinCohesion",
+ Dictionary_Entry_Value_FromDouble( -1.0e+04f ) ) );
+ contextExt->subDomainMaxCohesion = Dictionary_Entry_Value_AsDouble(
+ Dictionary_GetDefault( context->dictionary, "subDomainMaxCohesion",
+ Dictionary_Entry_Value_FromDouble( 3.0e+04f ) ) );
+
contextExt->trackLevel = Dictionary_Entry_Value_AsDouble(
Dictionary_GetDefault( context->dictionary, "trackLevel",
Dictionary_Entry_Value_FromDouble( 0.0f ) ) );
Modified: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h 2009-03-30 15:39:22 UTC (rev 14547)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Context.h 2009-03-30 19:26:02 UTC (rev 14548)
@@ -57,11 +57,14 @@
double xSubDomainFraction;
double ySubDomainFraction;
double zSubDomainFraction;
+ double xSubDomainLeft;
double xTriggerPointFraction;
double yTriggerPointFraction;
double zTriggerPointFraction;
double weakPointCohesion;
double triggerPointCohesion;
+ double subDomainMinCohesion;
+ double subDomainMaxCohesion;
double trackLevel;
double startThreshold;
double stopThreshold;
Modified: long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.c 2009-03-30 15:39:22 UTC (rev 14547)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/CreateWeakPoints.c 2009-03-30 19:26:02 UTC (rev 14548)
@@ -66,6 +66,7 @@
#endif
#define DEBUG
+//#define DEBUG2
/*
@@ -125,32 +126,41 @@
*/
if(fractionWeakPoints>=0.0 && fractionWeakPoints<=1.0) {
unsigned int rngSeed = contextExt->rngSeed;
- double x_subdomainFraction = contextExt->xSubDomainFraction;
- double y_subdomainFraction = contextExt->ySubDomainFraction;
- double z_subdomainFraction = contextExt->zSubDomainFraction;
+ double x_subDomainFraction = contextExt->xSubDomainFraction;
+ double y_subDomainFraction = contextExt->ySubDomainFraction;
+ double z_subDomainFraction = contextExt->zSubDomainFraction;
+ double x_subDomainLeftPosition = contextExt->xSubDomainLeft;
+ double subDomainMinCohesion = contextExt->subDomainMinCohesion;
+ double subDomainMaxCohesion = contextExt->subDomainMaxCohesion;
unsigned int numberSubDomainPoints,numberWeakPoints;
unsigned int *shuffleIndexListPtr, *randomNumberListPtr;
- int subdomain_I_element_range,subdomain_J_element_range,subdomain_K_element_range;
+ int subDomain_I_element_range,subDomain_J_element_range,subDomain_K_element_range,subDomain_I_element_left_position;
/* Report HillSlope plugin variables picked up (?) from xml parameter file */
Journal_Printf( context->snacInfo, "\tRNG seed = %u\n", rngSeed );
Journal_Printf( context->snacInfo, "\tFraction of weak points = %g\n", fractionWeakPoints );
- Journal_Printf( context->snacInfo, "\tx subdomain fraction = %g\n", x_subdomainFraction );
- Journal_Printf( context->snacInfo, "\ty subdomain fraction = %g\n", y_subdomainFraction );
- Journal_Printf( context->snacInfo, "\tz subdomain fraction = %g\n", z_subdomainFraction );
+ Journal_Printf( context->snacInfo, "\tx subDomain fraction = %g\n", x_subDomainFraction );
+ Journal_Printf( context->snacInfo, "\ty subDomain fraction = %g\n", y_subDomainFraction );
+ Journal_Printf( context->snacInfo, "\tz subDomain fraction = %g\n", z_subDomainFraction );
+ Journal_Printf( context->snacInfo, "\tx subDomain left position = %g\n", x_subDomainLeftPosition );
/*
* Define portion of mesh (cells) that need weak point "seeding"
* - lots of sloppy float/int casting back and forth here
*/
- subdomain_I_element_range = (int)(full_I_element_range*x_subdomainFraction);
- subdomain_J_element_range = (int)(full_J_element_range*y_subdomainFraction);
- subdomain_K_element_range = (int)(full_K_element_range*z_subdomainFraction);
+ subDomain_I_element_range = (int)(full_I_element_range*x_subDomainFraction);
+ subDomain_J_element_range = (int)(full_J_element_range*y_subDomainFraction);
+ subDomain_K_element_range = (int)(full_K_element_range*z_subDomainFraction);
+ if ( x_subDomainLeftPosition < 0.0 || subDomain_I_element_range+(int)(full_I_element_range*x_subDomainLeftPosition+0.5) > full_I_element_range)
+ subDomain_I_element_left_position = (full_I_element_range-subDomain_I_element_range)/2;
+ else
+ subDomain_I_element_left_position = (int)(full_I_element_range*x_subDomainLeftPosition+0.5);
- Journal_Printf( context->snacInfo, "\tx subdomain element range = %d/%d\n", subdomain_I_element_range,full_I_element_range );
- Journal_Printf( context->snacInfo, "\ty subdomain element range = %d/%d\n", subdomain_J_element_range,full_J_element_range );
- Journal_Printf( context->snacInfo, "\tz subdomain element range = %d/%d\n\n", subdomain_K_element_range,full_K_element_range );
+ Journal_Printf( context->snacInfo, "\tx subDomain element range = %d/%d\n", subDomain_I_element_range,full_I_element_range );
+ Journal_Printf( context->snacInfo, "\ty subDomain element range = %d/%d\n", subDomain_J_element_range,full_J_element_range );
+ Journal_Printf( context->snacInfo, "\tz subDomain element range = %d/%d\n", subDomain_K_element_range,full_K_element_range );
+ Journal_Printf( context->snacInfo, "\tx subDomain element left position = %d/%d\n\n", subDomain_I_element_left_position,full_I_element_range );
- numberSubDomainPoints = subdomain_I_element_range*subdomain_J_element_range*subdomain_K_element_range;
+ numberSubDomainPoints = subDomain_I_element_range*subDomain_J_element_range*subDomain_K_element_range;
numberWeakPoints = (unsigned int)((float)numberSubDomainPoints*fractionWeakPoints);
Journal_Printf(context->snacInfo, "Number of weak points = %u/%u\n", numberWeakPoints,numberSubDomainPoints);
@@ -164,16 +174,48 @@
* Create (1) ordered sequence of element indices, (2) parallel list of random variates
*/
index=0;
- for(index_I = 0; index_I < subdomain_I_element_range; index_I++) {
- for(index_J = 0; index_J < subdomain_J_element_range; index_J++) {
- for(index_K = 0; index_K < subdomain_K_element_range; index_K++) {
+ for(index_K = 0; index_K < subDomain_K_element_range; index_K++) {
+ for(index_J = 0; index_J < subDomain_J_element_range; index_J++) {
+ for(index_I = 0; index_I < subDomain_I_element_range; index_I++) {
/*
- * Work out the element index from the i,j,k and the required location of the subdomain
+ * Work out the element index from the i,j,k and the required location of the subDomain
*/
- element_gI = index_I+(full_I_element_range-subdomain_I_element_range)/2
+ element_gI = index_I + subDomain_I_element_left_position
+ full_I_element_range*(full_J_element_range-1-index_J)
- + full_I_element_range*full_J_element_range*(index_K+(full_K_element_range-subdomain_K_element_range)/2);
+ + full_I_element_range*full_J_element_range*(index_K+(full_K_element_range-subDomain_K_element_range)/2);
element_lI = Mesh_ElementMapGlobalToLocal( mesh, element_gI );
+
+
+ /*
+ * Impose heterogeneous cohesion within whole subDomain - if min/max values are non-neg
+ */
+ if(element_lI < mesh->elementLocalCount && subDomainMinCohesion>0.0) {
+ Snac_Element *element;
+ Snac_Material *material;
+ SnacPlastic_Element *plasticElement;
+ Tetrahedra_Index tetra_I;
+ double localRandomCohesion
+ = subDomainMinCohesion+((drand48()+drand48())/2.0)*(subDomainMaxCohesion-subDomainMinCohesion);
+ /*
+ * Naughty use of drand48 to generate triangular pdf variates spanning [subMin, subMax)
+ */
+
+ element = Snac_Element_At( context, element_lI );
+ material = &context->materialProperty[element->material_I];
+ plasticElement = ExtensionManager_Get( mesh->elementExtensionMgr, element, SnacPlastic_ElementHandle );
+ for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
+ plasticElement->plasticStrain[tetra_I] = PlasticStrainFromCohesion(material,localRandomCohesion);
+#ifdef DEBUG2
+ if(tetra_I==0)
+ Journal_Printf( context->snacInfo,
+ "Heterogeneous cohesion: timeStep=%d element (%d,%d,%d)#=%d->%d cohesion=%g plasticStrain=%g\n",
+ context->timeStep, index_I,index_J,index_K, element_gI, element_lI,
+ localRandomCohesion, plasticElement->plasticStrain[tetra_I] );
+#endif
+
+ }
+ } // End if
+
shuffleIndexListPtr[index]=element_lI;
randomNumberListPtr[index]=rand();
@@ -182,28 +224,30 @@
}
if(index>0 && index-1>numberSubDomainPoints) break;
}
- if(index>0 && index-1>numberSubDomainPoints) {fprintf(stderr, "We miscalculated the subdomain index\n"); break;}
+ if(index>0 && index-1>numberSubDomainPoints) {fprintf(stderr, "We miscalculated the subDomain index\n"); break;}
}
/*
* Sort the list of sub-domain local (ie global set pared down for this CPU)
* indices according to the random numbers
* - i.e., randomize the index list
*/
- ShellSort( randomNumberListPtr, shuffleIndexListPtr, numberSubDomainPoints );
+ if (numberWeakPoints>0) {
+ ShellSort( randomNumberListPtr, shuffleIndexListPtr, numberSubDomainPoints );
#ifdef DEBUG2
- for(index = 0; index < numberSubDomainPoints; index++) {
- Journal_Printf( context->snacInfo, "\t Subdomain point #%d = %d, RNG var=%d\n",
- index, shuffleIndexListPtr[index],randomNumberListPtr[index] );
+ for(index = 0; index < numberSubDomainPoints; index++) {
+ Journal_Printf( context->snacInfo, "\t SubDomain point #%d = %d, RNG var=%d\n",
+ index, shuffleIndexListPtr[index],randomNumberListPtr[index] );
+ }
+#endif
}
-#endif
/*
- * Now the subdomain index list is shuffled and the random number list is ordered
+ * Now the subDomain index list is shuffled and the random number list is ordered
*/
/*
- * Loop over the shuffled subdomain index list BUT only up to the number of weak points within it
- * - i.e., subset the subdomain and allocate weak points to only a random fraction of them
+ * Loop over the shuffled subDomain index list BUT only up to the number of weak points within it
+ * - i.e., subset the subDomain and allocate weak points to only a random fraction of them
*/
for(index = 0; index < numberWeakPoints; index++) {
element_lI = shuffleIndexListPtr[index];
@@ -227,7 +271,10 @@
for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
plasticElement->plasticStrain[tetra_I] = PlasticStrainFromCohesion(material,(double)contextExt->weakPointCohesion);
#ifdef DEBUG2
- if(tetra_I==0) Journal_Printf( context->snacInfo,"timeStep=%d et=%d %d plasticE: %e -> %e\n",context->timeStep, element_lI, tetra_I, material->plstrain[1], plasticElement->plasticStrain[tetra_I] );
+ if(tetra_I==0)
+ Journal_Printf( context->snacInfo,"Weak points: timeStep=%d element (%d,%d,%d)#=%d->%d cohesion=%g plasticStrain=%g\n",
+ context->timeStep, index_I,index_J,index_K, element_gI, element_lI,
+ contextExt->weakPointCohesion, plasticElement->plasticStrain[tetra_I] );
#endif
}
Modified: long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.c 2009-03-30 15:39:22 UTC (rev 14547)
+++ long/3D/SNAC/trunk/Snac/plugins/hillSlope/Register.c 2009-03-30 19:26:02 UTC (rev 14548)
@@ -108,15 +108,15 @@
SnacHillSlope_Type );
EntryPoint_Append(
- Context_GetEntryPoint( context, Snac_EP_LoopNodesMomentum ),
- SnacHillSlope_Type,
- SnacHillSlope_Track,
+ Context_GetEntryPoint( context, Snac_EP_Constitutive ),
+ SnacHillSlope_Type,
+ SnacHillSlope_CreateWeakPoints,
SnacHillSlope_Type );
EntryPoint_Append(
Context_GetEntryPoint( context, Snac_EP_LoopNodesMomentum ),
SnacHillSlope_Type,
- SnacHillSlope_CreateWeakPoints,
+ SnacHillSlope_Track,
SnacHillSlope_Type );
EntryPoint_Append(
Modified: long/3D/SNAC/trunk/Snac/plugins/plastic/Constitutive.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/plastic/Constitutive.c 2009-03-30 15:39:22 UTC (rev 14547)
+++ long/3D/SNAC/trunk/Snac/plugins/plastic/Constitutive.c 2009-03-30 19:26:02 UTC (rev 14548)
@@ -238,9 +238,11 @@
}
depls += plasticElement->plasticStrain[tetra_I]*element->tetra[tetra_I].volume;
totalVolume += element->tetra[tetra_I].volume;
+/* if(element_lI==738) fprintf(stderr," ***** %d: depm=%g acc plastic strain=%g\n", tetra_I, depm, element->tetra[tetra_I].volume); */
}
/* volume-averaged accumulated plastic strain, aps */
plasticElement->aps = depls/totalVolume;
+/* if(element_lI==738) fprintf(stderr," ***** acc plastic strain=%g\n", plasticElement->aps); */
}
}
More information about the CIG-COMMITS
mailing list