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

[becker at geodynamics.org]

Author: becker
Date: 2009-05-26 08:03:36 -0700 (Tue, 26 May 2009)
New Revision: 15053

Modified:
   mc/3D/CitcomS/trunk/lib/Initial_temperature.c
   mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
Log:
Fixed typo in blob temperature assignment. 

Changed conversion of coordinates, now assignment works. 



Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c	2009-05-26 09:41:14 UTC (rev 15052)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c	2009-05-26 15:03:36 UTC (rev 15053)
@@ -37,6 +37,7 @@
 void parallel_process_termination();
 void temperatures_conform_bcs(struct All_variables *);
 double modified_plgndr_a(int, int, double);
+void rtp2xyzd(double,double,double,double *);
 
 #include "initial_temperature.h"
 static void debug_tic(struct All_variables *);
@@ -136,7 +137,7 @@
 
     switch(E->convection.tic_method){
     case 2:			/* blob */
-      if( ! input_float_vector("blob_centqer", 3, E->convection.blob_center, m)) {
+      if( ! input_float_vector("blob_center", 3, E->convection.blob_center, m)) {
 	assert( E->sphere.caps == 12 || E->sphere.caps == 1 );
 	if(E->sphere.caps == 12) { /* Full version: just quit here */
 	  fprintf(stderr,"Missing input parameter: 'blob_center'.\n");
@@ -536,44 +537,38 @@
     int i, j ,k , m, node;
     int nox, noy, noz;
 
-    double theta_center, fi_center, r_center;
+    double theta_center, fi_center, r_center,x_center[4],dx[4];
     double radius, amp;
 
-    double x_center, y_center, z_center;
-    double x, y, z, distance;
+    double distance;
 
     noy = E->lmesh.noy;
     nox = E->lmesh.nox;
     noz = E->lmesh.noz;
 
     theta_center = E->convection.blob_center[0];
-    fi_center = E->convection.blob_center[1];
-    r_center = E->convection.blob_center[2];
-    radius = E->convection.blob_radius;
-    amp = E->convection.blob_dT;
+    fi_center    = E->convection.blob_center[1];
+    r_center     = E->convection.blob_center[2];
+    radius       = E->convection.blob_radius;
+    amp          = E->convection.blob_dT;
     
-    fprintf(stderr,"center=(%e %e %e) radius=%e dT=%e\n",
-            theta_center, fi_center, r_center, radius, amp);
+    if(E->parallel.me == 0)
+      fprintf(stderr,"center=(%e %e %e) radius=%e dT=%e\n",
+	      theta_center, fi_center, r_center, radius, amp);
+    
+    rtp2xyzd(r_center, theta_center, fi_center, (x_center+1));
 
-    x_center = r_center * sin(fi_center) * cos(theta_center);
-    y_center = r_center * sin(fi_center) * sin(theta_center);
-    z_center = r_center * cos(fi_center);
-
     /* compute temperature field according to nodal coordinate */
     for(m=1; m<=E->sphere.caps_per_proc; m++)
         for(i=1; i<=noy; i++)
             for(j=1; j<=nox;j ++)
                 for(k=1; k<=noz; k++) {
                     node = k + (j-1)*noz + (i-1)*nox*noz;
+		    dx[1] = E->x[m][1][node] - x_center[1];
+		    dx[2] = E->x[m][2][node] - x_center[2];
+		    dx[3] = E->x[m][3][node] - x_center[3];
+                    distance = sqrt(dx[1]*dx[1] + dx[2]*dx[2] + dx[3]*dx[3]);
 
-                    x = E->x[m][1][node];
-                    y = E->x[m][2][node];
-                    z = E->x[m][3][node];
-
-                    distance = sqrt((x-x_center)*(x-x_center) +
-                                    (y-y_center)*(y-y_center) +
-                                    (z-z_center)*(z-z_center));
-
                     if (distance < radius)
                         E->T[m][node] += amp * exp(-1.0*distance/radius);
                 }

Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c	2009-05-26 09:41:14 UTC (rev 15052)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c	2009-05-26 15:03:36 UTC (rev 15053)
@@ -53,6 +53,7 @@
 
 void calc_cbase_at_tp(float , float , float *);
 void rtp2xyz(float , float , float, float *);
+void rtp2xyzd(double , double , double, double *);
 void convert_pvec_to_cvec(float ,float , float , float *,float *);
 void *safe_malloc (size_t );
 void myerror(struct All_variables *,char *);
@@ -311,6 +312,15 @@
    sphere_to_cart
 
 */
+void rtp2xyzd(double r, double theta, double phi, double *xout)
+{
+  double rst;
+  rst = r * sin(theta);
+  xout[0] = rst * cos(phi);	/* x */
+  xout[1] = rst * sin(phi); 	/* y */
+  xout[2] = r * cos(theta);
+}
+/* float version */
 void rtp2xyz(float r, float theta, float phi, float *xout)
 {
   float rst;