[cig-commits] r13646 - in long/3D/Gale/trunk: . input/examples src/PICellerator/Utils/src
walter at geodynamics.org
walter at geodynamics.org
Thu Dec 11 01:29:58 PST 2008
Author: walter
Date: 2008-12-11 01:29:58 -0800 (Thu, 11 Dec 2008)
New Revision: 13646
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/input/examples/tibet.xml
long/3D/Gale/trunk/src/PICellerator/Utils/src/BuoyancyForceTerm.c
Log:
r2408 at dante: boo | 2008-12-11 01:21:29 -0800
Now apply density hydrostatic correction on a per-particle rather than per-element basis. This works for t>0 after mesh has been distorted. Currently hacked to only work with tibet model.
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2406
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2408
Modified: long/3D/Gale/trunk/input/examples/tibet.xml
===================================================================
--- long/3D/Gale/trunk/input/examples/tibet.xml 2008-12-11 09:15:35 UTC (rev 13645)
+++ long/3D/Gale/trunk/input/examples/tibet.xml 2008-12-11 09:29:58 UTC (rev 13646)
@@ -103,10 +103,6 @@
<param name="Rank">Scalar</param>
<param name="DataType">Double</param>
</struct>
- <struct name="pressureBCs">
- <param name="Type">CompositeVC</param>
- <param name="Data">mesh-linear</param>
- </struct>
<struct name="pressureDofLayout">
<param name="Type">DofLayout</param>
<param name="mesh">mesh-linear</param>
@@ -119,8 +115,6 @@
<param name="FEMesh">mesh-linear</param>
<param name="DofLayout">pressureDofLayout</param>
<param name="LinkedDofInfo">pressureLinkedDofs</param>
- <param name="BC">pressureBCs</param>
- <param name="IC">pressureICs</param>
</struct>
<struct name="backgroundDensity">
@@ -372,8 +366,6 @@
<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>
@@ -412,15 +404,6 @@
<param name="endZ">maxZ</param>
</struct>
-<!-- <struct name="stressBC"> -->
-<!-- <param name="Type">StressBC</param> -->
-<!-- <param name="ForceVector">mom_force</param> -->
-<!-- <param name="Swarm">picIntegrationPoints</param> -->
-<!-- <param name="wall">bottom</param> -->
-<!-- <param name="y_type">double</param> -->
-<!-- <param name="y_value">3237300000</param> -->
-<!-- </struct> -->
-
<struct name="surfaceAdaptor">
<param name="Type">SurfaceAdaptor</param>
<param name="mesh">mesh-linear</param>
@@ -459,6 +442,48 @@
<param name="MaterialPointsSwarm">materialSwarm</param>
</struct>
+
+<!-- <struct name="lowShape"> -->
+<!-- <param name="Type">Box</param> -->
+<!-- <param name="startX">450000</param> -->
+<!-- <param name="endX">550000</param> -->
+<!-- <param name="startY">68000</param> -->
+<!-- <param name="endY">80000</param> -->
+<!-- <param name="startZ">minZ</param> -->
+<!-- <param name="endZ">maxZ</param> -->
+<!-- </struct> -->
+<!-- <struct name="lowViscosity"> -->
+<!-- <param name="Type">FrankKamenetskii</param> -->
+<!-- <param name="TemperatureField">TemperatureField</param> -->
+ <!-- at T=273, viscosity=1e25, at T=1273, viscosity=1e20 -->
+<!-- <param name="eta0">2.31739465e-1</param> -->
+ <!-- <param name="eta0">2.31739465e26</param> -->
+<!-- <param name="theta">0.011512925465</param> -->
+<!-- </struct> -->
+<!-- <struct name="lowViscous"> -->
+<!-- <param name="Type">RheologyMaterial</param> -->
+<!-- <param name="Shape">lowShape</param> -->
+<!-- <param name="density">2800</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> -->
+<!-- <param name="Q">1.0e-12</param> -->
+<!-- <param name="lambda">0.0</param> -->
+<!-- </struct> -->
+<!-- </list> -->
+<!-- <list name="Rheology"> -->
+<!-- <param>lowViscosity</param> -->
+<!-- <param>yielding</param> -->
+<!-- <param>storeViscosity</param> -->
+<!-- <param>storeStress</param> -->
+<!-- </list> -->
+<!-- </struct> -->
+
+
+
<!-- Crust -->
<struct name="crustShape">
<param name="Type">PolygonShape</param>
@@ -511,10 +536,19 @@
<param name="frictionCoefficientAfterSoftening">0.57735026919</param>
</struct>
+<!-- <struct name="crustlowShape"> -->
+<!-- <param name="Type">Intersection</param> -->
+<!-- <list name="shapes"> -->
+<!-- <param>crustShape</param> -->
+<!-- <param>!lowShape</param> -->
+<!-- </list> -->
+<!-- </struct> -->
+
<struct name="crustViscous">
<param name="Type">RheologyMaterial</param>
+<!-- <param name="Shape">crustlowShape</param> -->
<param name="Shape">crustShape</param>
- <param name="density">3300</param>
+ <param name="density">2800</param>
<param name="alpha">3.0e-5</param>
<!-- <param name="density">upperDensity</param> -->
<!-- <param name="alpha">upperAlpha</param> -->
@@ -778,24 +812,6 @@
</list>
</struct>
- <struct name="pressureBCs">
- <param name="type">CompositeVC</param>
- <list name="vcList">
- <struct>
- <param name="type">WallVC</param>
- <param name="wall">top</param>
- <list name="variables">
- <struct>
- <param name="name">pressure</param>
- <param name="type">double</param>
- <param name="value">0</param>
- </struct>
- </list>
- </struct>
- </list>
- </struct>
-
-
<struct name="temperatureBCs">
<param name="type">CompositeVC</param>
<list name="vcList">
@@ -932,8 +948,7 @@
<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="upperDensity">2800</param>
<param name="lowerDensity">3300</param>
<param name="materialBoundary">68000</param>
<param name="gravity">9.81</param>
Modified: long/3D/Gale/trunk/src/PICellerator/Utils/src/BuoyancyForceTerm.c
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/BuoyancyForceTerm.c 2008-12-11 09:15:35 UTC (rev 13645)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/BuoyancyForceTerm.c 2008-12-11 09:29:58 UTC (rev 13646)
@@ -58,6 +58,45 @@
#include <assert.h>
#include <string.h>
+double DensityFunction(Coord coord)
+{
+ double T_0=273;
+ double A=0.0096667;
+ double B=93.333;
+ double C=1.0e-4;
+ double y_max=100000;
+ double material_boundary=68000;
+
+ double h=y_max-coord[1];
+ double alpha=3.0e-5;
+ double density, T;
+
+ if(h<0)
+ {
+ T=T_0;
+ }
+ else
+ {
+ T=T_0 + A*h + B*(1-exp(-C*h));
+ }
+ if(coord[1]>y_max)
+ {
+ density=0;
+ }
+ else if(coord[1]>material_boundary)
+ {
+ density=2800;
+ }
+ else
+ {
+ density=3300;
+ }
+
+ return density*(1-alpha*T);
+}
+
+
+
/* Textual name of this class */
const Type BuoyancyForceTerm_Type = "BuoyancyForceTerm";
@@ -346,96 +385,115 @@
void _BuoyancyForceTerm_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) {
- BuoyancyForceTerm* self = (BuoyancyForceTerm*) forceTerm;
- IntegrationPoint* particle;
- BuoyancyForceTerm_MaterialExt* materialExt;
- Particle_InCellIndex cParticle_I;
- Particle_InCellIndex cellParticleCount;
- Element_NodeIndex elementNodeCount;
- Dimension_Index dim = forceVector->dim;
- IntegrationPointsSwarm* swarm = (IntegrationPointsSwarm*)self->integrationSwarm;
- FeMesh* mesh = forceVector->feVariable->feMesh;
- Node_ElementLocalIndex eNode_I;
- Cell_Index cell_I;
- ElementType* elementType;
- Dof_Index nodeDofCount;
- double gravity;
- double detJac = 0.0;
- double factor;
- double Ni[8];
- double force;
- double* xi;
- Material* material;
- FeVariable* temperatureField = self->temperatureField;
- double temperature = 0.0;
- double *densityField = self->densityField;
- double background_density = 0.0;
- double* gHat;
- unsigned d_i;
+ BuoyancyForceTerm* self = (BuoyancyForceTerm*) forceTerm;
+ IntegrationPoint* particle;
+ BuoyancyForceTerm_MaterialExt* materialExt;
+ Particle_InCellIndex cParticle_I;
+ Particle_InCellIndex cellParticleCount;
+ Element_NodeIndex elementNodeCount;
+ Dimension_Index dim = forceVector->dim;
+ IntegrationPointsSwarm* swarm = (IntegrationPointsSwarm*)self->integrationSwarm;
+ FeMesh* mesh = forceVector->feVariable->feMesh;
+ Node_ElementLocalIndex eNode_I;
+ Cell_Index cell_I;
+ ElementType* elementType;
+ Dof_Index nodeDofCount;
+ double gravity;
+ double detJac = 0.0;
+ double factor;
+ double Ni[8];
+ double force;
+ double* xi;
+ Material* material;
+ FeVariable* temperatureField = self->temperatureField;
+ double temperature = 0.0;
+ double *densityField = self->densityField;
+ double background_density = 0.0;
+ double* gHat;
+ unsigned d_i;
+
+ double totalWeight = 0.0;
+ double adjustFactor = 0.0;
+
+ elementType = FeMesh_GetElementType( mesh, lElement_I );
+ elementNodeCount = elementType->nodeCount;
+ nodeDofCount = dim;
+ cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
+ cellParticleCount = swarm->cellParticleCountTbl[cell_I];
+ gHat = self->gHat;
+
+ /* adjust & adjustFactor -- 20060411 Alan
+ *
+ * The adjust decides whether an adjustFactor should be applied to
+ * the resulting factor. If on, the total weight of the
+ * particles in the cell are scaled to the cell local volume.
+ *
+ * This is designed to be used when integrating with swarms which do
+ * not cover the whole domain (ie - I used it to do dave.m's test of
+ * 1 swarm for blob, 1 swarm for background)
+ */
+ if ( self->adjust ) {
+ for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
+ particle = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, cell_I, cParticle_I );
+ totalWeight += particle->weight;
+ }
+ adjustFactor = swarm->weights->cellLocalVolume / totalWeight;
+ }
+ else {
+ adjustFactor = 1.0;
+ }
+
+ for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
+ Coord coord;
+
+ particle = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, cell_I, cParticle_I );
+ xi = particle->xi;
+
+ detJac = ElementType_JacobianDeterminant( elementType, mesh, lElement_I, xi, dim );
+ ElementType_EvaluateShapeFunctionsAt( elementType, xi, Ni );
+
+ /* Get parameters */
+ if ( temperatureField )
+ FeVariable_InterpolateWithinElement( temperatureField, lElement_I, xi, &temperature );
+/* if(densityField) */
+/* FeVariable_InterpolateWithinElement( densityField, lElement_I, xi, &background_density ); */
+
+ FeMesh_CoordLocalToGlobal(mesh, cell_I, xi, coord);
+ background_density=DensityFunction(coord);
+
+ 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=-gravity*(materialExt->density*(1.0-materialExt->alpha*temperature)
+ - background_density);
- double totalWeight = 0.0;
- double adjustFactor = 0.0;
-
- elementType = FeMesh_GetElementType( mesh, lElement_I );
- elementNodeCount = elementType->nodeCount;
- nodeDofCount = dim;
- cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
- cellParticleCount = swarm->cellParticleCountTbl[cell_I];
- gHat = self->gHat;
-
- /* adjust & adjustFactor -- 20060411 Alan
- *
- * The adjust decides whether an adjustFactor should be applied to the resulting factor.
- * If on, the total weight of the particles in the cell are scaled to the cell local volume.
- *
- * This is designed to be used when integrating with swarms which do not cover the whole domain
- * (ie - I used it to do dave.m's test of 1 swarm for blob, 1 swarm for background)
- */
- if ( self->adjust ) {
- for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
- particle = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, cell_I, cParticle_I );
- totalWeight += particle->weight;
- }
- adjustFactor = swarm->weights->cellLocalVolume / totalWeight;
- }
- else {
- adjustFactor = 1.0;
- }
-
- for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
- particle = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, cell_I, cParticle_I );
- xi = particle->xi;
-
- detJac = ElementType_JacobianDeterminant( elementType, mesh, lElement_I, xi, dim );
- ElementType_EvaluateShapeFunctionsAt( elementType, xi, Ni );
-
- /* 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 = - gravity*(materialExt->density * (1.0 - materialExt->alpha * temperature)
- - background_density);
-
- factor = detJac * particle->weight * adjustFactor * force;
-
- /* Apply force in verticle direction */
- for( eNode_I = 0 ; eNode_I < elementNodeCount; eNode_I++ ) {
- if( gHat ) {
- for( d_i = 0; d_i < dim; d_i++ )
- elForceVec[ eNode_I * nodeDofCount + d_i ] += gHat[d_i] * factor * Ni[ eNode_I ] ;
- }
- else {
- elForceVec[ eNode_I * nodeDofCount + J_AXIS ] += factor * Ni[ eNode_I ] ;
- }
- }
- }
-
+/* printf("force %d %d %g %g %g %g %g\n",lElement_I,cParticle_I,force,gravity, */
+/* materialExt->density*(1.0-materialExt->alpha*temperature), */
+/* background_density, */
+/* gravity*(materialExt->density*(1.0-materialExt->alpha*temperature) */
+/* - background_density)); */
+
+ factor = detJac * particle->weight * adjustFactor * force;
+
+ /* Apply force in verticle direction */
+ for( eNode_I = 0 ; eNode_I < elementNodeCount; eNode_I++ ) {
+ if( gHat ) {
+ for( d_i = 0; d_i < dim; d_i++ )
+ elForceVec[ eNode_I * nodeDofCount + d_i ] +=
+ gHat[d_i] * factor * Ni[ eNode_I ] ;
+ }
+ else {
+ elForceVec[ eNode_I * nodeDofCount + J_AXIS ] +=
+ factor * Ni[ eNode_I ] ;
+ }
+ }
+ }
}
/* The default implementation is for the gravity to be constant. */
More information about the CIG-COMMITS
mailing list