[cig-commits] r6583 - mc/3D/CitcomS/branches/compressible/lib

[tan2 at geodynamics.org]

Author: tan2
Date: 2007-04-16 15:37:27 -0700 (Mon, 16 Apr 2007)
New Revision: 6583

Modified:
   mc/3D/CitcomS/branches/compressible/lib/Advection_diffusion.c
   mc/3D/CitcomS/branches/compressible/lib/Size_does_matter.c
   mc/3D/CitcomS/branches/compressible/lib/global_defs.h
Log:
Compute Tmass matrix with density (missing heat capacity)

Modified: mc/3D/CitcomS/branches/compressible/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Advection_diffusion.c	2007-04-16 22:36:30 UTC (rev 6582)
+++ mc/3D/CitcomS/branches/compressible/lib/Advection_diffusion.c	2007-04-16 22:37:27 UTC (rev 6583)
@@ -347,7 +347,7 @@
     for (m=1;m<=E->sphere.caps_per_proc;m++)
         for(i=1;i<=E->lmesh.nno;i++) {
             if(!(E->node[m][i] & (TBX | TBY | TBZ)))
-                DTdot[m][i] *= E->Mass[m][i];         /* lumped mass matrix */
+                DTdot[m][i] *= E->TMass[m][i];         /* lumped mass matrix */
             else
                 DTdot[m][i] = 0.0;         /* lumped mass matrix */
         }

Modified: mc/3D/CitcomS/branches/compressible/lib/Size_does_matter.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Size_does_matter.c	2007-04-16 22:36:30 UTC (rev 6582)
+++ mc/3D/CitcomS/branches/compressible/lib/Size_does_matter.c	2007-04-16 22:37:27 UTC (rev 6583)
@@ -870,7 +870,7 @@
 void mass_matrix(struct All_variables *E)
 {
     int m,node,i,nint,e,lev;
-    int n[9];
+    int n[9], nz;
     void get_global_shape_fn();
     double myatan(),rtf[4][9],area,centre[4],temp[9],dx1,dx2,dx3;
     struct Shape_function GN;
@@ -959,8 +959,6 @@
                             *E->N.vpt[GNVINDEX(node,nint)];       /* int Na dV */
                 }
 
-                /* lumped mass matrix, equivalent to tmass in ConMan */
-                /* missing density*(heat capcity)? */
                 for(node=1;node<=enodes[E->mesh.nsd];node++)
                     E->MASS[lev][m][E->IEN[lev][m][e].node[node]] += temp[node];
 
@@ -983,6 +981,39 @@
     } /* end of for lev */
 
 
+    for (m=1;m<=E->sphere.caps_per_proc;m++) {
+
+        for(node=1;node<=E->lmesh.nno;node++)
+            E->TMass[m][node] = 0.0;
+
+        for(e=1;e<=E->lmesh.nel;e++)  {
+            get_global_shape_fn(E,e,&GN,&GNx,&dOmega,0,
+                                sphere_key,rtf,E->mesh.levmax,m);
+
+            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++) {
+                    /* XXX: missing heat capcity */
+                    temp[node] += E->refstate.rho[nz]
+                        * dOmega.vpt[nint]
+                        * g_point[nint].weight[E->mesh.nsd-1]
+                        * E->N.vpt[GNVINDEX(node,nint)]; /* int Na dV */
+                }
+            }
+
+            /* 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];
+
+        } /* end of for e */
+    } /* end of for m */
+
+    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];
+
+
     /* compute volume of this processor mesh and the whole mesh */
     E->lmesh.volume = 0;
     E->mesh.volume = 0;
@@ -1009,6 +1040,12 @@
                     fprintf(E->fp_out,"Mass[%d]= %g \n",node,E->MASS[lev][m][node]);
             }
         }
+
+        for (m=1;m<=E->sphere.caps_per_proc;m++)   {
+            fprintf(E->fp_out,"m=%d %d \n",E->sphere.capid[m],m);
+            for (node=1;node<=E->lmesh.nno;node++)
+                fprintf(E->fp_out,"TMass[%d]= %g \n",node,E->TMass[m][node]);
+        }
         fflush(E->fp_out);
     }
 

Modified: mc/3D/CitcomS/branches/compressible/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/global_defs.h	2007-04-16 22:36:30 UTC (rev 6582)
+++ mc/3D/CitcomS/branches/compressible/lib/global_defs.h	2007-04-16 22:37:27 UTC (rev 6583)
@@ -686,6 +686,7 @@
     double *temp[NCS],*temp1[NCS];
     float *NP[NCS],*edot[NCS],*Mass[NCS];
     float *MASS[MAX_LEVELS][NCS];
+    double *TMass[NCS];
     double *SX[MAX_LEVELS][NCS][4],*X[MAX_LEVELS][NCS][4];
     double *sx[NCS][4],*x[NCS][4];
     double *surf_det[NCS][5];