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

[tan2 at geodynamics.org]

Author: tan2
Date: 2007-09-07 15:53:22 -0700 (Fri, 07 Sep 2007)
New Revision: 7940

Modified:
   mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
   mc/3D/CitcomS/trunk/lib/Size_does_matter.c
Log:
Fixed tsolver.

Fixed a bug in TMass and another in convert element # to nz.


Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c	2007-09-07 01:38:52 UTC (rev 7939)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c	2007-09-07 22:53:22 UTC (rev 7940)
@@ -246,25 +246,31 @@
 	corrector(E,E->T,E->Tdot,DTdot);
 	temperatures_conform_bcs(E);
       }
-    }
 
-    if(E->advection.monitor_max_T) {
-        /* get the max temperature for new T */
-        E->monitor.T_interior = Tmaxd(E,E->T);
+      if(E->advection.monitor_max_T) {
+          /* get the max temperature for new T */
+          E->monitor.T_interior = Tmaxd(E,E->T);
 
-        /* if the max temperature changes too much, restore the old
-         * temperature field, calling the temperature solver using
-         * half of the timestep size */
-        if (E->monitor.T_interior/T_interior1 > E->monitor.T_maxvaried) {
-            for(m=1;m<=E->sphere.caps_per_proc;m++)
-                for (i=1;i<=E->lmesh.nno;i++)   {
-                    E->T[m][i] = T1[m][i];
-                    E->Tdot[m][i] = Tdot1[m][i];
-                }
-            iredo = 1;
-            E->advection.dt_reduced *= 0.5;
-            E->advection.last_sub_iterations ++;
-        }
+          /* if the max temperature changes too much, restore the old
+           * temperature field, calling the temperature solver using
+           * half of the timestep size */
+          if (E->monitor.T_interior/T_interior1 > E->monitor.T_maxvaried) {
+              if(E->parallel.me==0) {
+                  fprintf(stderr, "max T varied from %e to %e\n",
+                          T_interior1, E->monitor.T_interior);
+                  fprintf(E->fp, "max T varied from %e to %e\n",
+                          T_interior1, E->monitor.T_interior);
+              }
+              for(m=1;m<=E->sphere.caps_per_proc;m++)
+                  for (i=1;i<=E->lmesh.nno;i++)   {
+                      E->T[m][i] = T1[m][i];
+                      E->Tdot[m][i] = Tdot1[m][i];
+                  }
+              iredo = 1;
+              E->advection.dt_reduced *= 0.5;
+              E->advection.last_sub_iterations ++;
+          }
+      }
     }
 
   }  while ( iredo==1 && E->advection.last_sub_iterations <= 5);
@@ -601,7 +607,7 @@
       Q += E->composition.comp_el[m][0][el] * E->control.Q0ER;
     }
 
-    nz = ((el-1) % E->lmesh.noz) + 1;
+    nz = ((el-1) % E->lmesh.elz) + 1;
     rho = 0.5 * (E->refstate.rho[nz] + E->refstate.rho[nz+1]);
 
     if(E->control.disptn_number == 0)

Modified: mc/3D/CitcomS/trunk/lib/Size_does_matter.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Size_does_matter.c	2007-09-07 01:38:52 UTC (rev 7939)
+++ mc/3D/CitcomS/trunk/lib/Size_does_matter.c	2007-09-07 22:53:22 UTC (rev 7940)
@@ -990,27 +990,17 @@
             get_global_shape_fn(E,e,&GN,&GNx,&dOmega,0,
                                 sphere_key,rtf,E->mesh.levmax,m);
 
-            for(nint=1;nint<=vpts;nint++)
-                temp2[nint] = 0.0;
-
             for(node=1;node<=enodes[E->mesh.nsd];node++) {
+                temp[node] = 0.0;
                 nz = ((E->ien[m][e].node[node]-1) % E->lmesh.noz) + 1;
-                for(nint=1;nint<=vpts;nint++) {
-                    temp2[nint] = E->refstate.rho[nz]
+                for(nint=1;nint<=vpts;nint++)
+                    temp[node] += E->refstate.rho[nz]
                         * E->refstate.heat_capacity[nz]
+                        * dOmega.vpt[nint]
+                        * g_point[nint].weight[E->mesh.nsd-1]
                         * E->N.vpt[GNVINDEX(node,nint)];
-                }
             }
 
-            for(node=1;node<=enodes[E->mesh.nsd];node++) {
-                temp[node] = 0.0;
-                for(nint=1;nint<=vpts;nint++) {
-                    temp[node] += temp2[nint]
-                        * dOmega.vpt[nint]
-                        * g_point[nint].weight[E->mesh.nsd-1];
-                }
-            }
-
             /* lumped mass matrix, equivalent to tmass in ConMan */
             for(node=1;node<=enodes[E->mesh.nsd];node++)
                 E->TMass[m][E->ien[m][e].node[node]] += temp[node];
@@ -1018,6 +1008,7 @@
         } /* end of for e */
     } /* end of for m */
 
+    (E->exchange_node_d)(E,E->TMass,E->mesh.levmax);
     for (m=1;m<=E->sphere.caps_per_proc;m++)
         for(node=1;node<=E->lmesh.nno;node++)
             E->TMass[m][node] = 1.0 / E->TMass[m][node];