[cig-commits] r13640 - in long/3D/Gale/trunk: . input/cookbook input/examples src/PICellerator/Utils/src src/StgFEM/plugins/StandardConditionFunctions src/Underworld/Rheology/src
walter at geodynamics.org
walter at geodynamics.org
Wed Dec 10 04:32:15 PST 2008
Author: walter
Date: 2008-12-10 04:32:15 -0800 (Wed, 10 Dec 2008)
New Revision: 13640
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/input/cookbook/viscous_extension.xml
long/3D/Gale/trunk/input/examples/tibet.xml
long/3D/Gale/trunk/src/PICellerator/Utils/src/BuoyancyForceTerm.c
long/3D/Gale/trunk/src/PICellerator/Utils/src/BuoyancyForceTerm.h
long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.c
long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.h
long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h
Log:
r2404 at dante: boo | 2008-12-10 04:30:20 -0800
Added BackgroundDensityField and BackgroundPressureField to
subtract out the hydrostatic pressure that was causing
MixedStabiliser to go awry. Seems to work for viscous_extension
and tibet for a single step. I still need to make it work over
time.
Also, there is some debugging junk in the friction implementation
that I did not have time to get rid of.
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2399
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2404
Modified: long/3D/Gale/trunk/input/cookbook/viscous_extension.xml
===================================================================
--- long/3D/Gale/trunk/input/cookbook/viscous_extension.xml 2008-12-10 10:25:43 UTC (rev 13639)
+++ long/3D/Gale/trunk/input/cookbook/viscous_extension.xml 2008-12-10 12:32:15 UTC (rev 13640)
@@ -116,6 +116,59 @@
<param name="DofLayout">pressureDofLayout</param>
<param name="LinkedDofInfo">pressureLinkedDofs</param>
</struct>
+
+ <struct name="backgroundDensity">
+ <param name="Type">MeshVariable</param>
+ <param name="mesh">mesh-linear</param>
+ <param name="Rank">Scalar</param>
+ <param name="DataType">Double</param>
+ </struct>
+ <struct name="backgroundDensityICs">
+ <param name="Type">CompositeVC</param>
+ <param name="Data">mesh-linear</param>
+ </struct>
+ <struct name="backgroundDensityDofLayout">
+ <param name="Type">DofLayout</param>
+ <param name="mesh">mesh-linear</param>
+ <list name="BaseVariables">
+ <param>backgroundDensity</param>
+ </list>
+ </struct>
+ <struct name="backgroundDensityField">
+ <param name="Type">FeVariable</param>
+ <param name="FEMesh">mesh-linear</param>
+ <param name="DofLayout">backgroundDensityDofLayout</param>
+ <param name="LinkedDofInfo">backgroundDensityLinkedDofs</param>
+ <param name="BC">backgroundDensityBCs</param>
+ <param name="IC">backgroundDensityICs</param>
+ </struct>
+
+ <struct name="backgroundPressure">
+ <param name="Type">MeshVariable</param>
+ <param name="mesh">mesh-linear</param>
+ <param name="Rank">Scalar</param>
+ <param name="DataType">Double</param>
+ </struct>
+ <struct name="backgroundPressureICs">
+ <param name="Type">CompositeVC</param>
+ <param name="Data">mesh-linear</param>
+ </struct>
+ <struct name="backgroundPressureDofLayout">
+ <param name="Type">DofLayout</param>
+ <param name="mesh">mesh-linear</param>
+ <list name="BaseVariables">
+ <param>backgroundPressure</param>
+ </list>
+ </struct>
+ <struct name="BackgroundPressureField">
+ <param name="Type">FeVariable</param>
+ <param name="FEMesh">mesh-linear</param>
+ <param name="DofLayout">backgroundPressureDofLayout</param>
+ <param name="LinkedDofInfo">backgroundPressureLinkedDofs</param>
+ <param name="BC">backgroundPressureBCs</param>
+ <param name="IC">backgroundPressureICs</param>
+ </struct>
+
<struct name="StressField">
<param name="Type">StressField</param>
<param name="StrainRateField">StrainRateField</param>
@@ -317,12 +370,12 @@
<param name="storeVisc">storeViscosity</param>
<param name="StiffnessMatrix">c_matrix</param>
</struct>
-
<struct name="buoyancyForceTerm">
<param name="Type">BuoyancyForceTerm</param>
<param name="ForceVector">mom_force</param>
<param name="Swarm">picIntegrationPoints</param>
<param name="gravity">1.0</param>
+ <param name="BackgroundDensityField">backgroundDensityField</param>
</struct>
<struct name="background">
<param name="Type">Everywhere</param>
@@ -383,7 +436,7 @@
<param name="maxTimeSteps">10</param>
<param name="outputEvery">1</param>
<param name="dumpEvery">1</param>
- <param name="outputPath">./output.template</param>
+ <param name="outputPath">./output</param>
<param name="dim">2</param>
<param name="shadowDepth">1</param>
<param name="minX">0.0f</param>
@@ -485,22 +538,45 @@
</list>
</struct>
-
- <struct name="stressICs">
+ <struct name="backgroundDensityICs">
<param name="type">CompositeVC</param>
<list name="vcList">
<struct>
<param name="type">AllNodesVC</param>
<list name="variables">
<struct>
- <param name="name">stress</param>
+ <param name="name">backgroundDensity</param>
<param name="type">double</param>
- <param name="value">0.0</param>
+ <param name="value">1.0</param>
</struct>
</list>
- </struct>
+ </struct>
</list>
</struct>
+ <struct name="backgroundPressureICs">
+ <param name="type">CompositeVC</param>
+ <list name="vcList">
+ <struct>
+ <param name="type">AllNodesVC</param>
+ <list name="variables">
+ <struct>
+ <param name="name">backgroundPressure</param>
+ <param name="type">func</param>
+ <param name="value">Velocity_SimpleShear</param>
+ </struct>
+ </list>
+ </struct>
+ </list>
+ </struct>
+
+
<param name="checkpointEvery">1</param>
+ <param name="SimpleShearFactor">-1</param>
+ <param name="SimpleShearCentreY">0.35</param>
+
+ <param name="journal.info">True</param>
+ <param name="journal.debug">True</param>
+ <param name="journal-level.info">2</param>
+ <param name="journal-level.debug">2</param>
</StGermainData>
Modified: long/3D/Gale/trunk/input/examples/tibet.xml
===================================================================
--- long/3D/Gale/trunk/input/examples/tibet.xml 2008-12-10 10:25:43 UTC (rev 13639)
+++ long/3D/Gale/trunk/input/examples/tibet.xml 2008-12-10 12:32:15 UTC (rev 13640)
@@ -122,6 +122,57 @@
<param name="BC">pressureBCs</param>
<param name="IC">pressureICs</param>
</struct>
+
+ <struct name="backgroundDensity">
+ <param name="Type">MeshVariable</param>
+ <param name="mesh">mesh-linear</param>
+ <param name="Rank">Scalar</param>
+ <param name="DataType">Double</param>
+ </struct>
+ <struct name="backgroundDensityICs">
+ <param name="Type">CompositeVC</param>
+ <param name="Data">mesh-linear</param>
+ </struct>
+ <struct name="backgroundDensityDofLayout">
+ <param name="Type">DofLayout</param>
+ <param name="mesh">mesh-linear</param>
+ <list name="BaseVariables">
+ <param>backgroundDensity</param>
+ </list>
+ </struct>
+ <struct name="BackgroundDensityField">
+ <param name="Type">FeVariable</param>
+ <param name="FEMesh">mesh-linear</param>
+ <param name="DofLayout">backgroundDensityDofLayout</param>
+ <param name="LinkedDofInfo">backgroundDensityLinkedDofs</param>
+ <param name="IC">backgroundDensityICs</param>
+ </struct>
+
+ <struct name="backgroundPressure">
+ <param name="Type">MeshVariable</param>
+ <param name="mesh">mesh-linear</param>
+ <param name="Rank">Scalar</param>
+ <param name="DataType">Double</param>
+ </struct>
+ <struct name="backgroundPressureICs">
+ <param name="Type">CompositeVC</param>
+ <param name="Data">mesh-linear</param>
+ </struct>
+ <struct name="backgroundPressureDofLayout">
+ <param name="Type">DofLayout</param>
+ <param name="mesh">mesh-linear</param>
+ <list name="BaseVariables">
+ <param>backgroundPressure</param>
+ </list>
+ </struct>
+ <struct name="BackgroundPressureField">
+ <param name="Type">FeVariable</param>
+ <param name="FEMesh">mesh-linear</param>
+ <param name="DofLayout">backgroundPressureDofLayout</param>
+ <param name="LinkedDofInfo">backgroundPressureLinkedDofs</param>
+ <param name="IC">backgroundPressureICs</param>
+ </struct>
+
<struct name="StressField">
<param name="Type">StressField</param>
<param name="StrainRateField">StrainRateField</param>
@@ -321,24 +372,18 @@
<param name="Swarm">gaussSwarm</param>
<param name="picSwarm">picIntegrationPoints</param>
<param name="storeVisc">storeViscosity</param>
+ <param name="height">maxY</param>
+<!-- <param name="scale_factor">stabiliserScale</param> -->
<param name="StiffnessMatrix">c_matrix</param>
</struct>
-<!-- <struct name="hydrostaticCorrection"> -->
-<!-- <param name="Type">HydrostaticCorrection</param> -->
-<!-- <param name="ForceVector">cont_force</param> -->
-<!-- <param name="Swarm">gaussSwarm</param> -->
-<!-- <param name="picSwarm">picIntegrationPoints</param> -->
-<!-- <param name="storeVisc">storeViscosity</param> -->
-<!-- <param name="force_type">func</param> -->
-<!-- <param name="force_value">PressureProfile</param> -->
-<!-- </struct> -->
<struct name="buoyancyForceTerm">
<param name="Type">BuoyancyForceTerm</param>
<param name="ForceVector">mom_force</param>
-<!-- <param name="TemperatureField">TemperatureField</param> -->
+ <param name="TemperatureField">TemperatureField</param>
+ <param name="BackgroundDensityField">BackgroundDensityField</param>
<param name="Swarm">picIntegrationPoints</param>
- <param name="gravity">9.81</param>
+ <param name="gravity">gravity</param>
</struct>
<struct name="background">
<param name="Type">Everywhere</param>
@@ -455,6 +500,7 @@
<struct name="yielding">
<param name="Type">MohrCoulomb</param>
<param name="PressureField">PressureField</param>
+ <param name="BackgroundPressureField">BackgroundPressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
<param name="Context">context</param>
@@ -468,9 +514,10 @@
<struct name="crustViscous">
<param name="Type">RheologyMaterial</param>
<param name="Shape">crustShape</param>
-<!-- <param name="density">3300</param> -->
- <param name="density">2800</param>
+ <param name="density">3300</param>
<param name="alpha">3.0e-5</param>
+<!-- <param name="density">upperDensity</param> -->
+<!-- <param name="alpha">upperAlpha</param> -->
<param name="diffusivity">1.0e-6</param>
<list name="heatingElements">
<struct>
@@ -508,6 +555,8 @@
<param name="Shape">mantleShape</param>
<param name="density">3300</param>
<param name="alpha">3.0e-5</param>
+<!-- <param name="density">lowerDensity</param> -->
+<!-- <param name="alpha">lowerAlpha</param> -->
<param name="diffusivity">1.0e-6</param>
<list name="heatingElements">
<struct>
@@ -837,6 +886,38 @@
</list>
</struct>
+ <struct name="backgroundDensityICs">
+ <param name="type">CompositeVC</param>
+ <list name="vcList">
+ <struct>
+ <param name="type">AllNodesVC</param>
+ <list name="variables">
+ <struct>
+ <param name="name">backgroundDensity</param>
+ <param name="type">func</param>
+ <param name="value">DensityProfile</param>
+ </struct>
+ </list>
+ </struct>
+ </list>
+ </struct>
+
+ <struct name="backgroundPressureICs">
+ <param name="type">CompositeVC</param>
+ <list name="vcList">
+ <struct>
+ <param name="type">AllNodesVC</param>
+ <list name="variables">
+ <struct>
+ <param name="name">backgroundPressure</param>
+ <param name="type">func</param>
+ <param name="value">PressureProfile</param>
+ </struct>
+ </list>
+ </struct>
+ </list>
+ </struct>
+
<!-- These numbers have to be kept in sync with the top and bottom BC's. The
parser can not read the variables and use them there, so they have
to be specified in two places. -->
@@ -849,13 +930,17 @@
<!-- 1/h_r = 10^-4 -->
<param name="TemperatureProfileExponentialCoefficient2">1.0e-4</param>
+ <param name="upperAlpha">3.0e-5</param>
+ <param name="lowerAlpha">3.0e-5</param>
+ <param name="upperDensity">3300</param>
+<!-- <param name="upperDensity">2800</param> -->
+ <param name="lowerDensity">3300</param>
+ <param name="materialBoundary">68000</param>
+ <param name="gravity">9.81</param>
+
<!-- <param name="checkpointEvery">1</param> -->
<param name="dtFactor">1.0</param>
-<!-- <param name="SimpleShearCentreY">100000</param> -->
-<!-- <param name="SimpleShearFactor">-31547.8105236</param> -->
- <param name="SimpleShearCentreY">101510.391518</param>
- <param name="SimpleShearFactor">-31078.4049316</param>
<param name="journal.info">True</param>
<param name="journal.debug">True</param>
<param name="journal-level.info">2</param>
Modified: long/3D/Gale/trunk/src/PICellerator/Utils/src/BuoyancyForceTerm.c
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/BuoyancyForceTerm.c 2008-12-10 10:25:43 UTC (rev 13639)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/BuoyancyForceTerm.c 2008-12-10 12:32:15 UTC (rev 13640)
@@ -66,6 +66,7 @@
ForceVector* forceVector,
Swarm* integrationSwarm,
FeVariable* temperatureField,
+ FeVariable* densityField,
double gravity,
Bool adjust,
Materials_Register* materials_Register )
@@ -77,6 +78,7 @@
forceVector,
integrationSwarm,
temperatureField,
+ densityField,
gravity,
adjust,
materials_Register );
@@ -129,11 +131,13 @@
void _BuoyancyForceTerm_Init(
BuoyancyForceTerm* self,
FeVariable* temperatureField,
+ FeVariable* densityField,
double gravity,
Bool adjust,
Materials_Register* materials_Register )
{
self->temperatureField = temperatureField;
+ self->densityField = densityField;
self->gravity = gravity;
self->gHat = NULL;
self->adjust = adjust;
@@ -145,6 +149,7 @@
ForceVector* forceVector,
Swarm* integrationSwarm,
FeVariable* temperatureField,
+ FeVariable* densityField,
double gravity,
Bool adjust,
Materials_Register* materials_Register )
@@ -152,7 +157,8 @@
BuoyancyForceTerm* self = (BuoyancyForceTerm*) forceTerm;
ForceTerm_InitAll( self, forceVector, integrationSwarm, NULL );
- _BuoyancyForceTerm_Init( self, temperatureField, gravity, adjust, materials_Register );
+ _BuoyancyForceTerm_Init( self, temperatureField, densityField,
+ gravity, adjust, materials_Register );
}
void _BuoyancyForceTerm_Delete( void* forceTerm ) {
@@ -177,6 +183,7 @@
/* General info */
Journal_PrintPointer( stream, self->temperatureField );
+ Journal_PrintPointer( stream, self->densityField );
Journal_PrintDouble( stream, self->gravity );
}
@@ -203,7 +210,7 @@
Dictionary* dict;
Dictionary_Entry_Value* tmp;
char* rootKey;
- FeVariable* temperatureField;
+ FeVariable *temperatureField, *densityField;
double gravity;
Bool adjust;
Materials_Register* materials_Register;
@@ -217,6 +224,8 @@
dict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( cf->componentDict, self->name ) );
temperatureField = Stg_ComponentFactory_ConstructByKey( cf, self->name, "TemperatureField", FeVariable, False, data ) ;
+ densityField = Stg_ComponentFactory_ConstructByKey( cf, self->name, "BackgroundDensityField",
+ FeVariable, False, data ) ;
gravity = Stg_ComponentFactory_GetDouble( cf, self->name, "gravity", 0.0 );
adjust = Stg_ComponentFactory_GetBool( cf, self->name, "adjust", False );
@@ -243,7 +252,8 @@
materials_Register = Stg_ObjectList_Get( cf->registerRegister, "Materials_Register" );
assert( materials_Register );
- _BuoyancyForceTerm_Init( self, temperatureField, gravity, adjust, materials_Register );
+ _BuoyancyForceTerm_Init( self, temperatureField, densityField,
+ gravity, adjust, materials_Register );
}
void _BuoyancyForceTerm_Build( void* forceTerm, void* data ) {
@@ -261,6 +271,8 @@
if ( self->temperatureField )
Stg_Component_Build( self->temperatureField, data, False );
+ if ( self->densityField )
+ Stg_Component_Build( self->densityField, data, False );
/* Sort out material extension stuff */
self->materialExtHandle = Materials_Register_AddMaterialExtension(
@@ -315,6 +327,8 @@
if ( self->temperatureField )
Stg_Component_Initialise( self->temperatureField, data, False );
+ if ( self->densityField )
+ Stg_Component_Initialise( self->densityField, data, False );
for ( i = 0; i < self->materialSwarmCount; ++i ) {
Stg_Component_Initialise( self->densitySwarmVariables[i], data, False );
@@ -354,6 +368,8 @@
Material* material;
FeVariable* temperatureField = self->temperatureField;
double temperature = 0.0;
+ double *densityField = self->densityField;
+ double background_density = 0.0;
double* gHat;
unsigned d_i;
@@ -396,12 +412,16 @@
/* Get parameters */
if ( temperatureField )
FeVariable_InterpolateWithinElement( temperatureField, lElement_I, xi, &temperature );
+ if(densityField)
+ FeVariable_InterpolateWithinElement( densityField, lElement_I, xi, &background_density );
material = IntegrationPointsSwarm_GetMaterialOn( (IntegrationPointsSwarm*) swarm, particle );
materialExt = ExtensionManager_Get( material->extensionMgr, material, self->materialExtHandle );
/* Calculate Force */
gravity = BuoyancyForceTerm_CalcGravity( self, (Swarm*)swarm, lElement_I, particle );
- force = - materialExt->density * gravity * (1.0 - materialExt->alpha * temperature);
+ force = - gravity*(materialExt->density * (1.0 - materialExt->alpha * temperature)
+ - background_density);
+
factor = detJac * particle->weight * adjustFactor * force;
/* Apply force in verticle direction */
Modified: long/3D/Gale/trunk/src/PICellerator/Utils/src/BuoyancyForceTerm.h
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/BuoyancyForceTerm.h 2008-12-10 10:25:43 UTC (rev 13639)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/BuoyancyForceTerm.h 2008-12-10 12:32:15 UTC (rev 13640)
@@ -68,6 +68,7 @@
\
/* BuoyancyForceTerm info */ \
FeVariable* temperatureField; \
+ FeVariable* densityField; \
double gravity; \
double* gHat; \
Bool adjust; \
@@ -84,6 +85,7 @@
ForceVector* forceVector,
Swarm* integrationSwarm,
FeVariable* temperatureField,
+ FeVariable* densityField,
double gravity,
Bool adjust,
Materials_Register* materials_Register );
@@ -109,6 +111,7 @@
ForceVector* forceVector,
Swarm* integrationSwarm,
FeVariable* temperatureField,
+ FeVariable* densityField,
double gravity,
Bool adjust,
Materials_Register* materials_Register );
Modified: long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.c 2008-12-10 10:25:43 UTC (rev 13639)
+++ long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.c 2008-12-10 12:32:15 UTC (rev 13640)
@@ -131,6 +131,12 @@
condFunc = ConditionFunction_New( StgFEM_StandardConditionFunctions_TemperatureProfile, "TemperatureProfile");
ConditionFunction_Register_Add( context->condFunc_Register, condFunc );
+ condFunc = ConditionFunction_New( StgFEM_StandardConditionFunctions_PressureProfile, "PressureProfile");
+ ConditionFunction_Register_Add( context->condFunc_Register, condFunc );
+
+ condFunc = ConditionFunction_New( StgFEM_StandardConditionFunctions_DensityProfile, "DensityProfile");
+ ConditionFunction_Register_Add( context->condFunc_Register, condFunc );
+
condFunc = ConditionFunction_New( StG_FEM_StandardConditionFunctions_Gaussian, "Gaussian");
ConditionFunction_Register_Add( context->condFunc_Register, condFunc );
@@ -1466,6 +1472,113 @@
}
}
+double PressureProfile_Analytic(double density, double gravity, double h,
+ double alpha, double T_0, double A,
+ double B, double C)
+{
+ return density*gravity*h*(1-alpha*(T_0 + A*h/2 + B*(1+exp(-C*h)/(C*h))));
+}
+
+void StgFEM_StandardConditionFunctions_PressureProfile( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) {
+ FiniteElementContext * context = (FiniteElementContext*)_context;
+ FeVariable* feVariable = NULL;
+ FeMesh* mesh = NULL;
+ Dictionary* dictionary = context->dictionary;
+ double* result = (double*) _result;
+ double* coord;
+ double T_0, A, B, C;
+ double upper_density, lower_density, upper_alpha, lower_alpha, y_max,
+ material_boundary;
+ double gravity;
+
+ feVariable = (FeVariable*)FieldVariable_Register_GetByName( context->fieldVariable_Register, "VelocityField" );
+ mesh = feVariable->feMesh;
+ coord = Mesh_GetVertex( mesh, node_lI );
+
+ T_0 = Dictionary_GetDouble( dictionary, "TemperatureProfileTop");
+ A = Dictionary_GetDouble( dictionary, "TemperatureProfileLinearCoefficient");
+ B = Dictionary_GetDouble( dictionary,
+ "TemperatureProfileExponentialCoefficient1");
+ C = Dictionary_GetDouble( dictionary,
+ "TemperatureProfileExponentialCoefficient2");
+ y_max = Dictionary_GetDouble( dictionary, "maxY");
+ upper_density=Dictionary_GetDouble(dictionary,"upperDensity");
+ lower_density=Dictionary_GetDouble(dictionary,"lowerDensity");
+ upper_alpha=Dictionary_GetDouble( dictionary, "upperAlpha");
+ lower_alpha=Dictionary_GetDouble( dictionary, "lowerAlpha");
+ material_boundary=Dictionary_GetDouble(dictionary, "materialBoundary");
+ gravity=Dictionary_GetDouble( dictionary, "gravity" );
+ if(coord[1]>=y_max)
+ {
+ *result=0;
+ }
+ else
+ {
+ double h=y_max-coord[1];
+ if(coord[1]>material_boundary)
+ {
+ *result=PressureProfile_Analytic(upper_density, gravity, h,
+ upper_alpha, T_0, A, B, C);
+ }
+ else
+ {
+ *result=
+ PressureProfile_Analytic(upper_density, gravity, material_boundary,
+ upper_alpha, T_0, A, B, C)
+ - PressureProfile_Analytic(lower_density, gravity,material_boundary,
+ lower_alpha, T_0, A, B, C)
+ + PressureProfile_Analytic(lower_density, gravity,h,
+ lower_alpha, T_0, A, B, C);
+ }
+
+/* printf("pressure %d %g %g %g %g %g\n",node_lI,coord[1],h, */
+/* material_boundary, */
+
+ }
+}
+
+void StgFEM_StandardConditionFunctions_DensityProfile( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) {
+ FiniteElementContext * context = (FiniteElementContext*)_context;
+ FeVariable* feVariable = NULL;
+ FeMesh* mesh = NULL;
+ Dictionary* dictionary = context->dictionary;
+ double* result = (double*) _result;
+ double* coord;
+ double upper_density, lower_density, upper_alpha, lower_alpha, y_max,
+ material_boundary;
+
+ feVariable = (FeVariable*)FieldVariable_Register_GetByName( context->fieldVariable_Register, "VelocityField" );
+ mesh = feVariable->feMesh;
+ coord = Mesh_GetVertex( mesh, node_lI );
+
+ y_max = Dictionary_GetDouble_WithDefault( dictionary, "maxY", 0.0 );
+ upper_density=Dictionary_GetDouble_WithDefault(dictionary,"upperDensity",0.0);
+
+ lower_density=Dictionary_GetDouble_WithDefault(dictionary,"lowerDensity",0.0);
+ upper_alpha=Dictionary_GetDouble_WithDefault( dictionary, "upperAlpha", 0.0 );
+ lower_alpha=Dictionary_GetDouble_WithDefault( dictionary, "lowerAlpha", 0.0 );
+ material_boundary=Dictionary_GetDouble_WithDefault(dictionary,
+ "materialBoundary",0.0);
+ if(coord[1]>y_max)
+ {
+ *result=0;
+ }
+ else
+ {
+ double T;
+ StgFEM_StandardConditionFunctions_TemperatureProfile(node_lI,var_I,
+ _context,&T);
+ if(coord[1]>material_boundary)
+ {
+ *result=upper_density*(1-upper_alpha*T);
+ }
+ else
+ {
+ *result=lower_density*(1-lower_alpha*T);
+ }
+ }
+}
+
void StgFEM_StandardConditionFunctions_ERF( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) {
FiniteElementContext * context = (FiniteElementContext*)_context;
FeVariable* feVariable = NULL;
@@ -1566,7 +1679,7 @@
stressField = (FeVariable*)FieldVariable_Register_GetByName
( context->fieldVariable_Register, "StressField" );
pressureField = (FeVariable*)FieldVariable_Register_GetByName
- ( context->fieldVariable_Register, "NodalPressureField" );
+ ( context->fieldVariable_Register, "PressureField" );
viscosityField = (FeVariable*)FieldVariable_Register_GetByName
( context->fieldVariable_Register, "ViscosityField" );
@@ -1613,73 +1726,127 @@
{
*result=boundary_initial_v;
+ if(inds[tangent_dim]<8 && tangent_dim==0)
+ *result=0;
+
printf("friction %d %d %g initial\n",tangent_dim,inds[tangent_dim],*result);
}
- else if(*(context->nonLinearIteration_I)==0
- || (fabs(tangent_stress)<yield_stress && v[tangent_dim]==boundary_v))
+ else if(context->nonLinearIteration_I!=NULL
+ && (*(context->nonLinearIteration_I)!=0
+ || context->timeStep>1)
+ && fabs(tangent_stress)<yield_stress
+ && v[tangent_dim]==boundary_v)
{
*result=boundary_v;
/* *result=v[tangent_dim]; */
- printf("friction %d %d %g %g old\n",tangent_dim,inds[tangent_dim],
- *result,tangent_stress);
+ printf("friction %d %d %g %g %g %g %g old\n",
+ tangent_dim,inds[tangent_dim],
+ *result,v[0],v[1],boundary_v,tangent_stress);
}
else
{
- Node_LocalIndex temp_node;
- double *temp_coord, delta, v_off[3];
+/* FeVariable* dvField = NULL; */
+/* double dv[9]; */
+/* dvField = (FeVariable*)FieldVariable_Register_GetByName */
+/* ( context->fieldVariable_Register, "VelocityGradientsField" ); */
+/* FeVariable_GetValueAtNode(dvField,node_lI,dv); */
- if(inds[normal_dim]==0)
- ++inds[normal_dim];
- else
- --inds[normal_dim];
-
- node_gI=RegularMeshUtils_Node_3DTo1D(feMesh, inds);
- Mesh_GlobalToDomain(feMesh, MT_VERTEX, node_gI, &temp_node);
-
- FeVariable_GetValueAtNode(velocityField,temp_node,v_off);
- temp_coord=Mesh_GetVertex(feMesh, temp_node);
-
- delta=fabs(temp_coord[normal_dim]-coord[normal_dim]);
-
- /* Make sure that we always make the velocity closer to the
- boundary velocity */
- if(v_off[tangent_dim]>boundary_v)
- *result=v_off[tangent_dim]-delta*yield_stress/viscosity;
- else
- *result=v_off[tangent_dim]+delta*yield_stress/viscosity;
-
- /* Average the results so that we don't oscillate a lot, looking
- for a solution. */
- *result=(*result+v[tangent_dim])/2;
-
- /* If the yield stress is negative (e.g. the pressure is
- negative), then the boundary moves in lockstep with the
- element above it. */
- if(yield_stress<0)
+/* printf("friction %d %d %g %g %g %g dv\n",tangent_dim,inds[tangent_dim], */
+/* dv[0],dv[1],dv[2],dv[3]); */
+
+
+ /* If we have never yielded before, or if the tangent stress is
+ very different from the yield stress, then use the gradient to
+ compute a first approximation. */
+
+/* if(v[tangent_dim]==boundary_v */
+/* || fabs((yield_stress-fabs(tangent_stress)) */
+/* /(yield_stress+fabs(tangent_stress)))>0.5 */
+/* || yield_stress<0) */
{
- *result=v_off[tangent_dim];
+ Node_LocalIndex temp_node;
+ double *temp_coord, delta, v_off[3];
- printf("friction %d %d %g %g set\n",tangent_dim,inds[tangent_dim],
- *result, v_off[tangent_dim]);
- }
+ if(inds[normal_dim]==0)
+ ++inds[normal_dim];
+ else
+ --inds[normal_dim];
+
+ node_gI=RegularMeshUtils_Node_3DTo1D(feMesh, inds);
+ Mesh_GlobalToDomain(feMesh, MT_VERTEX, node_gI, &temp_node);
+
+ FeVariable_GetValueAtNode(velocityField,temp_node,v_off);
+ temp_coord=Mesh_GetVertex(feMesh, temp_node);
+
+ delta=fabs(temp_coord[normal_dim]-coord[normal_dim]);
+
+ /* Make sure that we always make the velocity closer to the
+ boundary velocity */
+ if(v[tangent_dim]>boundary_v)
+ *result=v[tangent_dim]
+ - delta*(yield_stress-fabs(tangent_stress))/viscosity;
+ else
+ *result=v[tangent_dim]
+ + delta*(yield_stress-fabs(tangent_stress))/viscosity;
- /* If we overshot, make the node stick */
- if((*result-boundary_v)*(v_off[tangent_dim]-boundary_v)<0)
- {
- printf("friction %d %d %g %g %g overshoot\n",
- tangent_dim,inds[tangent_dim],
- *result, boundary_v, v_off[tangent_dim]);
- *result=boundary_v;
+/* if(v_off[tangent_dim]>boundary_v) */
+/* *result=v_off[tangent_dim]-delta*yield_stress/viscosity; */
+/* else */
+/* *result=v_off[tangent_dim]+delta*yield_stress/viscosity; */
+
+ /* Average the results so that we don't oscillate a lot, looking
+ for a solution. */
+/* *result=(*result+v[tangent_dim])/2; */
+
+ /* If the yield stress is negative (e.g. the pressure is
+ negative), then the boundary moves in lockstep with the
+ element above it. */
+ if(yield_stress<0)
+ {
+ *result=v_off[tangent_dim];
+
+ printf("friction %d %d %g %g set\n",tangent_dim,inds[tangent_dim],
+ *result, v_off[tangent_dim]);
+ }
+
+ /* If we overshot, make the node stick */
+ if((*result-boundary_v)*(v_off[tangent_dim]-boundary_v)<0)
+ {
+ printf("friction %d %d %g %g %g overshoot\n",
+ tangent_dim,inds[tangent_dim],
+ *result, boundary_v, v_off[tangent_dim]);
+ *result=boundary_v + (v_off[tangent_dim]-boundary_v)*1.0e-5;
+/* *result=boundary_v; */
+ }
+
+
+ printf("friction %d %d %g %g %g %g %g %g %g %g %g initiate\n",
+ tangent_dim, inds[tangent_dim],
+ *result,delta,
+ yield_stress,normal_stress,p,v_off[tangent_dim],tangent_stress,
+ viscosity,v[tangent_dim]);
}
+/* else */
+/* { */
+/* *result=boundary_v + (v[tangent_dim]-boundary_v) */
+/* *(1+2*(yield_stress-fabs(tangent_stress))/(yield_stress+fabs(tangent_stress))); */
+
+/* /\* If we overshot, make the node stick *\/ */
+/* if((*result-boundary_v)*(v[tangent_dim]-boundary_v)<0) */
+/* { */
+/* printf("friction %d %d %g %g %g overshoot\n", */
+/* tangent_dim,inds[tangent_dim], */
+/* *result, boundary_v, v[tangent_dim]); */
+/* *result=boundary_v; */
+/* } */
-
- printf("friction %d %d %g %g %g %g %g %g %g slide\n",
- tangent_dim, inds[tangent_dim],
- *result,delta,
- yield_stress,normal_stress,p,v_off[tangent_dim],tangent_stress);
-
-
+/* printf("friction %d %d %g %g %g %g %g %g %g slide\n",tangent_dim, */
+/* inds[tangent_dim],*result, */
+/* yield_stress,tangent_stress,v[tangent_dim],boundary_v, */
+/* (yield_stress-tangent_stress),(yield_stress+tangent_stress) */
+/* ); */
+/* } */
}
}
Modified: long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.h 2008-12-10 10:25:43 UTC (rev 13639)
+++ long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.h 2008-12-10 12:32:15 UTC (rev 13640)
@@ -88,6 +88,8 @@
void StG_FEM_StandardConditionFunctions_StepFunctionProduct3( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;
void StG_FEM_StandardConditionFunctions_StepFunctionProduct4( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;
void StgFEM_StandardConditionFunctions_TemperatureProfile( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;
+void StgFEM_StandardConditionFunctions_PressureProfile( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;
+void StgFEM_StandardConditionFunctions_DensityProfile( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;
void StG_FEM_StandardConditionFunctions_Gaussian( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;
void StgFEM_StandardConditionFunctions_ConvectionBenchmark( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2008-12-10 10:25:43 UTC (rev 13639)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2008-12-10 12:32:15 UTC (rev 13640)
@@ -110,6 +110,7 @@
void _MohrCoulomb_Init(
MohrCoulomb* self,
FeVariable* pressureField,
+ FeVariable* backgroundPressureField,
FeVariable* velocityGradientsField,
MaterialPointsSwarm* materialPointsSwarm,
double cohesion,
@@ -120,6 +121,7 @@
{
self->materialPointsSwarm = materialPointsSwarm;
self->pressureField = pressureField;
+ self->backgroundPressureField = backgroundPressureField;
self->velocityGradientsField = velocityGradientsField;
self->cohesion = cohesion;
@@ -160,6 +162,7 @@
void *data ){
MohrCoulomb* self = (MohrCoulomb*)rheology;
FeVariable* pressureField;
+ FeVariable* backgroundPressureField;
MaterialPointsSwarm* materialPointsSwarm;
FeVariable* velocityGradientsField;
@@ -177,12 +180,15 @@
"MaterialPointsSwarm", MaterialPointsSwarm, True, data );
pressureField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
"PressureField", FeVariable, True, data );
+ backgroundPressureField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
+ "BackgroundPressureField", FeVariable, False, data );
velocityGradientsField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
"VelocityGradientsField", FeVariable, True, data );
_MohrCoulomb_Init(
self,
pressureField,
+ backgroundPressureField,
velocityGradientsField,
materialPointsSwarm,
Stg_ComponentFactory_GetDouble( cf, self->name, "cohesion", 0.0 ),
@@ -197,6 +203,8 @@
/* Build parent */
_YieldRheology_Build( self, data );
+ if ( self->backgroundPressureField )
+ Stg_Component_Build( self->backgroundPressureField, data, False );
}
void _MohrCoulomb_Initialise( void* rheology, void* data ) {
@@ -207,6 +215,9 @@
_YieldRheology_Initialise( self, data );
+ if ( self->backgroundPressureField )
+ Stg_Component_Initialise( self->backgroundPressureField, data, False );
+
/* Initialise these components now just in case this function is called before their own _Initialise function */
Stg_Component_Initialise( self->materialPointsSwarm, data, False );
Stg_Component_Initialise( self->strainWeakening, data, False );
@@ -296,12 +307,19 @@
double effectiveFrictionCoefficient;
double frictionalStrength;
- FeVariable* pressureField = self->pressureField;
- double pressure;
+ FeVariable* pressureField = self->pressureField;
+ FeVariable* backgroundPressureField = self->backgroundPressureField;
+ double pressure, backgroundPressure;
double theta;
double factor;
FeVariable_InterpolateWithinElement( pressureField, lElement_I, xi, &pressure );
+ if(backgroundPressureField)
+ {
+ FeVariable_InterpolateWithinElement( backgroundPressureField, lElement_I, xi,
+ &backgroundPressure );
+ pressure+=backgroundPressure;
+ }
/* Calculate frictional strength. We modify the friction and
cohesion because we have grouped terms from the normal
@@ -388,10 +406,15 @@
Dimension_Index dim = constitutiveMatrix->dim;
Eigenvector evectors[3];
double e0, e1, e2;
- double trace;
+ double trace, p;
int i;
FeVariable_InterpolateWithinElement( self->pressureField, lElement_I, xi, &self->currentPressure );
+ if(self->backgroundPressureField)
+ {
+ FeVariable_InterpolateWithinElement( self->backgroundPressureField, lElement_I, xi, &p );
+ self->currentPressure+=p;
+ }
FeVariable_InterpolateWithinElement( self->velocityGradientsField, lElement_I, xi, self->currentVelocityGradient );
TensorArray_GetSymmetricPart( self->currentVelocityGradient, dim, strainRate );
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h 2008-12-10 10:25:43 UTC (rev 13639)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h 2008-12-10 12:32:15 UTC (rev 13640)
@@ -59,6 +59,7 @@
/* Param passed in */ \
MaterialPointsSwarm* materialPointsSwarm; \
FeVariable* pressureField; \
+ FeVariable* backgroundPressureField; \
FeVariable* velocityGradientsField; \
/* Director component is used to update the normal */\
Director* director; \
More information about the CIG-COMMITS
mailing list