[cig-commits] [commit] rajesh-petsc-schur: Changed E->P to be 0-based; results dont match master (f0aeb23)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/citcoms

On branch  : rajesh-petsc-schur
Link       : https://github.com/geodynamics/citcoms/compare/464e1b32299b15819f93efd98d969cddb84dfe51...f97ae655a50bdbd6dac1923a3471ee4dae178fbd

>---------------------------------------------------------------

commit f0aeb23394aee331871b91cd102e3091d8f28576
Author: Rajesh Kommu <rajesh.kommu at gmail.com>
Date:   Wed Sep 10 11:36:26 2014 -0700

    Changed E->P to be 0-based; results dont match master


>---------------------------------------------------------------

f0aeb23394aee331871b91cd102e3091d8f28576
 lib/Checkpoints.c  |  4 ++--
 lib/Instructions.c |  6 +++---
 lib/Topo_gravity.c | 10 +++++-----
 3 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/lib/Checkpoints.c b/lib/Checkpoints.c
index 24b073e..d8618c8 100644
--- a/lib/Checkpoints.c
+++ b/lib/Checkpoints.c
@@ -460,7 +460,7 @@ static void momentum_checkpoint(struct All_variables *E, FILE *fp)
      * and won't be used when read it. */
     for(m=1; m<=E->sphere.caps_per_proc; m++) {
         /* Pressure at equation points */
-        fwrite(E->P[m], sizeof(double), E->lmesh.npno+1, fp);
+        fwrite(E->P[m], sizeof(double), E->lmesh.npno, fp);
 
         /* velocity at equation points */
         fwrite(E->U[m], sizeof(double), E->lmesh.neq, fp);
@@ -488,7 +488,7 @@ static void read_momentum_checkpoint(struct All_variables *E, FILE *fp)
 
     for(m=1; m<=E->sphere.caps_per_proc; m++) {
         /* Pressure at equation points */
-      if(fread(E->P[m], sizeof(double), E->lmesh.npno+1, fp) !=  E->lmesh.npno+1)
+      if(fread(E->P[m], sizeof(double), E->lmesh.npno, fp) !=  E->lmesh.npno+1)
 	myerror(E,"read_momentum_checkpoint: error at P");
         /* velocity at equation points */
       if(fread(E->U[m], sizeof(double), E->lmesh.neq, fp) != E->lmesh.neq)
diff --git a/lib/Instructions.c b/lib/Instructions.c
index 1e7b3d6..6269f76 100644
--- a/lib/Instructions.c
+++ b/lib/Instructions.c
@@ -953,7 +953,7 @@ void allocate_common_vars(E)
   elx  = E->lmesh.elx;
   ely  = E->lmesh.ely;
 
-  E->P[j]        = (double *) malloc((npno+1)*sizeof(double));
+  E->P[j]        = (double *) malloc(npno*sizeof(double));
   E->T[j]        = (double *) malloc((nno+1)*sizeof(double));
   E->NP[j]       = (float *) malloc((nno+1)*sizeof(float));
   E->buoyancy[j] = (double *) malloc((nno+1)*sizeof(double));
@@ -1149,7 +1149,7 @@ void allocate_common_vars(E)
       E->heating_latent[j][i] = 1.0;
   }
 
-  for(i=1;i<=E->lmesh.npno;i++)
+  for(i=0;i<E->lmesh.npno;i++)
       E->P[j][i] = 0.0;
 
   mat_prop_allocate(E);
@@ -1451,7 +1451,7 @@ void initial_pressure(E)
     report(E,"Initialize pressure field");
 
   for (m=1;m<=E->sphere.caps_per_proc;m++)
-    for(i=1;i<=E->lmesh.npno;i++)
+    for(i=0;i<E->lmesh.npno;i++)
       E->P[m][i]=0.0;
 
   return;
diff --git a/lib/Topo_gravity.c b/lib/Topo_gravity.c
index 2c1ea7e..5481db6 100644
--- a/lib/Topo_gravity.c
+++ b/lib/Topo_gravity.c
@@ -432,9 +432,9 @@ void compute_nodal_stress(struct All_variables *E,
       vor /= E->eco[m][e].area;
 
       /* add the pressure term */
-      Szz -= E->P[m][e];
-      Sxx -= E->P[m][e];
-      Syy -= E->P[m][e];
+      Szz -= E->P[m][e-1]; /* E->P is now 0-based. Change e-1 to e */
+      Sxx -= E->P[m][e-1]; /* when the loop variable e goes from */
+      Syy -= E->P[m][e-1]; /* 0 to nel-1, instead of the current 1 to nel */
 
       for(i=1;i<=ends;i++) {
         node = E->ien[m][e].node[i]; /* assign to global nodes */
@@ -1084,8 +1084,8 @@ void get_CBF_topo(E,H,HB)       /* call this only for top and bottom processors*
 //      get_elt_g(E,el,eltg,lev,j);
 
       for(m=0;m<dims*ends;m++) {
-           res[m]  = eltf[m]  - E->elt_del[lev][j][el].g[m][0]  * E->P[j][el];
-           resb[m] = eltfb[m] - E->elt_del[lev][j][elb].g[m][0]* E->P[j][elb];
+           res[m]  = eltf[m]  - E->elt_del[lev][j][el].g[m][0]  * E->P[j][el-1];
+           resb[m] = eltfb[m] - E->elt_del[lev][j][elb].g[m][0]* E->P[j][elb-1];
 //           res[m]  = eltf[m] - eltg[m][0]  * E->P[j][el];
 //           resb[m] = eltfb[m] - eltgb[m][0]* E->P[j][elb];
             }