[cig-commits] r14746 - in long/3D/Gale/trunk: . src/StgFEM/plugins/StandardConditionFunctions

[walter at geodynamics.org]

Author: walter
Date: 2009-04-17 15:20:02 -0700 (Fri, 17 Apr 2009)
New Revision: 14746

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.c
Log:
 r2654 at dante:  boo | 2009-04-16 06:38:49 -0700
 Added Garrett Ito's generalization of the temperature profile



Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
   - 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2652
   + 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2654

Modified: long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.c	2009-04-17 20:30:21 UTC (rev 14745)
+++ long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.c	2009-04-17 22:20:02 UTC (rev 14746)
@@ -1416,32 +1416,61 @@
 
 
 void StgFEM_StandardConditionFunctions_TemperatureProfile( 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, y_max;
-
-	feVariable = (FeVariable*)FieldVariable_Register_GetByName( context->fieldVariable_Register, "VelocityField" );
-	mesh       = feVariable->feMesh;
-	coord      = Mesh_GetVertex( mesh, node_lI );
-
-	T_0 = Dictionary_GetDouble_WithDefault( dictionary, "TemperatureProfileTop", 0.0 );
-	A = Dictionary_GetDouble_WithDefault( dictionary, "TemperatureProfileLinearCoefficient", 0.0 );
-	B = Dictionary_GetDouble_WithDefault( dictionary, "TemperatureProfileExponentialCoefficient1", 0.0 );
-	C = Dictionary_GetDouble_WithDefault( dictionary, "TemperatureProfileExponentialCoefficient2", 0.0 );
-        y_max = Dictionary_GetDouble_WithDefault( dictionary, "maxY", 0.0 );
-
-        if(coord[1]>y_max)
-          {
-            *result=T_0;
-          }
-        else
-          {
-            *result=T_0 + A*(y_max-coord[1]) + B*(1-exp(-C*(y_max-coord[1])));
-          }
+  FiniteElementContext *	context            = (FiniteElementContext*)_context;
+  FeVariable*             feVariable         = NULL;
+  FeMesh*     mesh               = NULL;
+  Dictionary*             dictionary         = context->dictionary;
+  double*                 result             = (double*) _result;
+  double*                 coord;
+  double                  T_0, H_0, dH, H, H_m, A, B, C, x_min, x_max, y_max, T_m, xc, dum;
+  /* G.Ito 10/08 added variables x_min, x_max, T_m, Xc, to do variation in x
+     and limit maximum T */
+  
+  feVariable = (FeVariable*)FieldVariable_Register_GetByName( context->fieldVariable_Register, "VelocityField" );
+  mesh       = feVariable->feMesh;
+  coord      = Mesh_GetVertex( mesh, node_lI );
+  
+  T_0 = Dictionary_GetDouble_WithDefault( dictionary, "TemperatureProfileTop", 0.0 );
+  T_m = Dictionary_GetDouble_WithDefault( dictionary, "TemperatureProfileMax", 10000.0 );
+  H_0 = Dictionary_GetDouble_WithDefault( dictionary, "TemperatureProfileH0", -1.0 );
+  H_m = Dictionary_GetDouble_WithDefault( dictionary, "TemperatureProfileHm", 1.0e+8 );
+  dH = Dictionary_GetDouble_WithDefault( dictionary, "TemperatureProfiledH", 0.0 );     
+  A = Dictionary_GetDouble_WithDefault( dictionary, "TemperatureProfileLinearCoefficient", 0.0 );
+  B = Dictionary_GetDouble_WithDefault( dictionary, "TemperatureProfileExponentialCoefficient1", 0.0 );
+  C = Dictionary_GetDouble_WithDefault( dictionary, "TemperatureProfileExponentialCoefficient2", 0.0 );
+  y_max = Dictionary_GetDouble_WithDefault( dictionary, "maxY", 0.0 );
+  x_max = Dictionary_GetDouble_WithDefault( dictionary, "maxX", 0.0 );
+  x_min = Dictionary_GetDouble_WithDefault( dictionary, "minX", 0.0 );
+  xc = Dictionary_GetDouble_WithDefault( dictionary, "ExtensionCentreX", 0.0 );
+  
+  if (H_0<0.0)
+    {
+      if(coord[1]>y_max)
+        {
+          *result=T_0;
+        }
+      else
+        {
+          *result=T_0 + A*(y_max-coord[1]) + B*(1-exp(-C*(y_max-coord[1])));
+        }
+    }
+  else
+    {
+      if(coord[1]>=y_max)
+        {
+          *result=T_0;
+        }
+      else
+        {
+          H=H_0 + 2*fabs(coord[0]-xc)/(x_max-x_min)*dH;
+          if (H>H_m) H=H_m;
+          
+          dum=T_0 + ((T_m-T_0)/H)*(y_max-coord[1])
+            + B*(1-exp(-C*(y_max-coord[1])));
+          if (dum>T_m) dum=T_m;
+          *result=dum;
+        }
+    }
 }
 
 void StgFEM_StandardConditionFunctions_ERF( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) {