[cig-commits] r15054 - mc/3D/CitcomS/trunk/lib

[becker at geodynamics.org]

Author: becker
Date: 2009-05-26 08:21:42 -0700 (Tue, 26 May 2009)
New Revision: 15054

Modified:
   mc/3D/CitcomS/trunk/lib/Initial_temperature.c
   mc/3D/CitcomS/trunk/lib/convection_variables.h
Log:
Added blob_bc_persist parameter (default = off) to allow assigning
blob-like temperature anomalies to BCs.



Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c	2009-05-26 15:03:36 UTC (rev 15053)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c	2009-05-26 15:21:42 UTC (rev 15054)
@@ -134,7 +134,7 @@
     input_float("half_space_age", &(E->convection.half_space_age), "40.0,1e-3,nomax", m);
     input_float("mantle_temp",&(E->control.mantle_temp),"1.0",m);
 
-
+    
     switch(E->convection.tic_method){
     case 2:			/* blob */
       if( ! input_float_vector("blob_center", 3, E->convection.blob_center, m)) {
@@ -152,6 +152,7 @@
       }
       input_float("blob_radius", &(E->convection.blob_radius), "0.063,0.0,1.0", m);
       input_float("blob_dT", &(E->convection.blob_dT), "0.18,nomin,nomax", m);
+      input_boolean("blob_bc_persist",&(E->convection.blob_bc_persist),"off",m);
       break;
     case 4:
       /*
@@ -538,14 +539,18 @@
     int nox, noy, noz;
 
     double theta_center, fi_center, r_center,x_center[4],dx[4];
-    double radius, amp;
-
+    double radius, amp, r1,rout,rin;
+    const double e_4 = 1e-4;
     double distance;
 
     noy = E->lmesh.noy;
     nox = E->lmesh.nox;
     noz = E->lmesh.noz;
 
+    rout = E->sphere.ro;
+    rin = E->sphere.ri;
+
+
     theta_center = E->convection.blob_center[0];
     fi_center    = E->convection.blob_center[1];
     r_center     = E->convection.blob_center[2];
@@ -569,8 +574,19 @@
 		    dx[3] = E->x[m][3][node] - x_center[3];
                     distance = sqrt(dx[1]*dx[1] + dx[2]*dx[2] + dx[3]*dx[3]);
 
-                    if (distance < radius)
-                        E->T[m][node] += amp * exp(-1.0*distance/radius);
+                    if (distance < radius){
+		      E->T[m][node] += amp * exp(-1.0*distance/radius);
+
+		      if(E->convection.blob_bc_persist){
+			r1 = E->sx[m][3][node];
+			if((fabs(r1 - rout) < e_4) || (fabs(r1 - rin) < e_4)){
+			  /* at bottom or top of box, assign as TBC */
+			  E->sphere.cap[m].TB[1][node]=E->T[m][node];
+			  E->sphere.cap[m].TB[2][node]=E->T[m][node];
+			  E->sphere.cap[m].TB[3][node]=E->T[m][node];
+			}
+		      }
+		    }
                 }
     return;
 }

Modified: mc/3D/CitcomS/trunk/lib/convection_variables.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/convection_variables.h	2009-05-26 15:03:36 UTC (rev 15053)
+++ mc/3D/CitcomS/trunk/lib/convection_variables.h	2009-05-26 15:21:42 UTC (rev 15054)
@@ -42,6 +42,7 @@
   float blob_center[3];
   float blob_radius;
   float blob_dT;
+  int blob_bc_persist;
 
     struct SOURCES {
 	    int number;