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

[becker at geodynamics.org]

Author: becker
Date: 2008-11-01 11:37:20 -0700 (Sat, 01 Nov 2008)
New Revision: 13213

Modified:
   mc/3D/CitcomS/trunk/lib/Full_sphere_related.c
   mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
   mc/3D/CitcomS/trunk/lib/Instructions.c
   mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
Log:
Fixed two typos, velocity norm function was called where pressure norm should have been called. 
This leads to core dumps, all seems to be working now. 



Modified: mc/3D/CitcomS/trunk/lib/Full_sphere_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_sphere_related.c	2008-11-01 18:13:44 UTC (rev 13212)
+++ mc/3D/CitcomS/trunk/lib/Full_sphere_related.c	2008-11-01 18:37:20 UTC (rev 13213)
@@ -202,6 +202,8 @@
   qx = malloc((temp+1)*sizeof(double));
   qy = malloc((temp+1)*sizeof(double));
 
+
+
   /* 4 corners of the cap in Cartesian coordinates */
   /* the order of corners is: */
   /*  1 - 4 */
@@ -242,6 +244,7 @@
                             center[2]*center[2]));;
   referencep[1] = myatan(center[1], center[0]);
 
+
   lev = E->mesh.levmax;
 
      /* # of elements/nodes per cap */
@@ -416,12 +419,12 @@
        lvnox = E->lmesh.NOX[lev];
        lvnoy = E->lmesh.NOY[lev];
        lvnoz = E->lmesh.NOZ[lev];
-       
+
        node = 1;
        for (k=0; k<lvnoy; k++) {
 	 for (j=0, ns=step*lnoy*k; j<lvnox; j++, ns+=step) {
 	   theta = qx[ns];
-	   fi = qy[ns];
+	   fi =    qy[ns];
 	   
 	   cost = cos(theta);
 	   sint = sin(theta);
@@ -435,9 +438,9 @@
 	     E->SX[lev][m][3][node] = E->sphere.R[lev][i];
 	     
 	     /*   x,y,and z oordinates   */
-	     E->X[lev][m][1][node] = E->sphere.R[lev][i]*sint*cosf;
-	     E->X[lev][m][2][node] = E->sphere.R[lev][i]*sint*sinf;
-	     E->X[lev][m][3][node] = E->sphere.R[lev][i]*cost;
+	     E->X[lev][m][1][node]  = E->sphere.R[lev][i]*sint*cosf;
+	     E->X[lev][m][2][node]  = E->sphere.R[lev][i]*sint*sinf;
+	     E->X[lev][m][3][node]  = E->sphere.R[lev][i]*cost;
 	     
 	     node++;
 	   }
@@ -474,6 +477,7 @@
   free ((void *)qx    );
   free ((void *)qy    );
 
+
   return;
 }
 

Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c	2008-11-01 18:13:44 UTC (rev 13212)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c	2008-11-01 18:37:20 UTC (rev 13213)
@@ -189,7 +189,6 @@
 
   void full_coord_of_cap();
   void compute_angle_surf_area ();
-
   rr = (double *)  malloc((E->mesh.noz+1)*sizeof(double));
   RR = (double *)  malloc((E->mesh.noz+1)*sizeof(double));
 
@@ -259,7 +258,6 @@
      full_coord_of_cap(E,j,ii);
      }
 
-
   if (E->control.verbose) {
       for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++)   {
           fprintf(E->fp_out,"output_coordinates before rotation %d \n",lev);

Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c	2008-11-01 18:13:44 UTC (rev 13212)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c	2008-11-01 18:37:20 UTC (rev 13213)
@@ -114,6 +114,7 @@
 
            /* physical domain */
     (E->solver.node_locations)(E);
+    if(chatty)fprintf(stderr,"node locations done\n");
 
     allocate_velocity_vars(E);
     if(chatty)fprintf(stderr,"velocity vars done\n");

Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c	2008-11-01 18:13:44 UTC (rev 13212)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c	2008-11-01 18:37:20 UTC (rev 13213)
@@ -286,7 +286,7 @@
                                    E->monitor.incompressibility);
     }
 
-
+  
     r0dotz0 = 0;
 
     while( (count < *steps_max) &&
@@ -304,7 +304,6 @@
         r1dotz1 = global_pdot(E, r1, z1, lev);
         assert(r1dotz1 != 0.0  /* Division by zero in head of incompressibility iteration */);
 
-
         /* update search direction */
         if(count == 0)
             for (m=1; m<=E->sphere.caps_per_proc; m++)
@@ -318,7 +317,6 @@
                     s2[m][j] = z1[m][j] + delta * s1[m][j];
         }
 
-
         /* solve K*u1 = grad(s2) for u1 */
         assemble_grad_p(E, s2, F, lev);
         valid = solve_del2_u(E, E->u1, F, imp*v_res, lev);
@@ -361,7 +359,7 @@
         v_norm = sqrt(E->monitor.vdotv);
         p_norm = sqrt(E->monitor.pdotp);
         dvelocity = alpha * sqrt(global_v_norm2(E, E->u1) / (1e-32 + E->monitor.vdotv));
-        dpressure = alpha * sqrt(global_v_norm2(E, s2) / (1e-32 + E->monitor.pdotp));
+        dpressure = alpha * sqrt(global_p_norm2(E, s2) / (1e-32 + E->monitor.pdotp));
 
         /* how many consecutive converging iterations? */
         if(dvelocity < imp && dpressure < imp)
@@ -647,7 +645,7 @@
         v_norm = sqrt(E->monitor.vdotv);
         p_norm = sqrt(E->monitor.pdotp);
         dvelocity = sqrt(global_v_norm2(E, F) / (1e-32 + E->monitor.vdotv));
-        dpressure = sqrt(global_v_norm2(E, s0) / (1e-32 + E->monitor.pdotp));
+        dpressure = sqrt(global_p_norm2(E, s0) / (1e-32 + E->monitor.pdotp));
 
         /* how many consecutive converging iterations? */
         if(dvelocity < imp && dpressure < imp)