[cig-commits] r11908 - mc/3D/CitcomCU/trunk/src

[tan2 at geodynamics.org]

Author: tan2
Date: 2008-05-02 16:56:00 -0700 (Fri, 02 May 2008)
New Revision: 11908

Modified:
   mc/3D/CitcomCU/trunk/src/Construct_arrays.c
   mc/3D/CitcomCU/trunk/src/Element_calculations.c
   mc/3D/CitcomCU/trunk/src/General_matrix_functions.c
   mc/3D/CitcomCU/trunk/src/Instructions.c
   mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c
Log:
Fixed some memory problems in the Stokes solver.

Some arrays were not allcated with enough elements (off by 1). Writing to the
end of those arrays can potentially overwrite other parts of the memory. In
practice, it seems that this bug didn't affect the solver at all, but this is
a subtle bug to detect.

This bug was in CitcomS pre-1.0 version and is likely in other versions of
citcom family as well.


Modified: mc/3D/CitcomCU/trunk/src/Construct_arrays.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Construct_arrays.c	2008-05-02 23:43:25 UTC (rev 11907)
+++ mc/3D/CitcomCU/trunk/src/Construct_arrays.c	2008-05-02 23:56:00 UTC (rev 11908)
@@ -906,7 +906,7 @@
 
 	for(level = E->mesh.levmax; level >= E->mesh.levmin; level--)
 	{
-		for(j = 0; j < E->lmesh.NEQ[level]; j++)
+		for(j = 0; j <= E->lmesh.NEQ[level] + 1; j++)
 			E->temp[j] = 0.0;
 
 		for(i = 1; i <= E->lmesh.NNO[level]; i++)

Modified: mc/3D/CitcomCU/trunk/src/Element_calculations.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Element_calculations.c	2008-05-02 23:43:25 UTC (rev 11907)
+++ mc/3D/CitcomCU/trunk/src/Element_calculations.c	2008-05-02 23:56:00 UTC (rev 11908)
@@ -377,9 +377,10 @@
 	//const int dofs = E->mesh.dof;
 	const int max_eqn = max_eqn_interaction[dims];
 
-	for(e = 0; e <= neq; e++)
+	for(e = 0; e <= neq + 1; e++)
 		Au[e] = 0.0;
 
+	u[neq + 1] = 0;
 	loc0 = 1;
 
 	for(e = 1; e <= nno; e++)

Modified: mc/3D/CitcomCU/trunk/src/General_matrix_functions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/General_matrix_functions.c	2008-05-02 23:43:25 UTC (rev 11907)
+++ mc/3D/CitcomCU/trunk/src/General_matrix_functions.c	2008-05-02 23:56:00 UTC (rev 11908)
@@ -684,9 +684,9 @@
 		r2 = (double *)malloc((1 + E->lmesh.NEQ[mem_lev]) * sizeof(double));
 		z0 = (double *)malloc((1 + E->lmesh.NEQ[mem_lev]) * sizeof(double));
 		z1 = (double *)malloc((1 + E->lmesh.NEQ[mem_lev]) * sizeof(double));
-		p1 = (double *)malloc((1 + E->lmesh.NEQ[mem_lev]) * sizeof(double));
-		p2 = (double *)malloc((1 + E->lmesh.NEQ[mem_lev]) * sizeof(double));
-		Ap = (double *)malloc((1 + E->lmesh.NEQ[mem_lev]) * sizeof(double));
+		p1 = (double *)malloc((2 + E->lmesh.NEQ[mem_lev]) * sizeof(double));
+		p2 = (double *)malloc((2 + E->lmesh.NEQ[mem_lev]) * sizeof(double));
+		Ap = (double *)malloc((2 + E->lmesh.NEQ[mem_lev]) * sizeof(double));
 		been_here++;
 	}
 
@@ -1103,13 +1103,13 @@
 		}
 
 
-	temp = (higher_precision *) malloc((neq + 1) * sizeof(higher_precision));
+	temp = (higher_precision *) malloc((neq + 2) * sizeof(higher_precision));
 
 	count = 0;
 
 	while(count < steps)
 	{
-		for(i = 0; i <= neq; i++)
+		for(i = 0; i <= neq + 1; i++)
 			temp[i] = zeroo;
 		for(i = 1; i <= nno; i++)
 		{
@@ -1196,14 +1196,14 @@
 
 	if(been_here == 0)
 	{
-		temp = (higher_precision *) malloc((E->lmesh.NEQ[E->mesh.levmax] + 1) * sizeof(higher_precision));
+		temp = (higher_precision *) malloc((E->lmesh.NEQ[E->mesh.levmax] + 2) * sizeof(higher_precision));
 		temp1 = (higher_precision *) malloc((E->lmesh.NEQ[E->mesh.levmax] + 1) * sizeof(higher_precision));
 		been_here++;
 	}
 
 	while(count < steps)
 	{
-		for(i = 0; i <= neq; i++)
+		for(i = 0; i <= neq + 1; i++)
 			temp[i] = zeroo;
 
 		if(dims == 3)

Modified: mc/3D/CitcomCU/trunk/src/Instructions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Instructions.c	2008-05-02 23:43:25 UTC (rev 11907)
+++ mc/3D/CitcomCU/trunk/src/Instructions.c	2008-05-02 23:56:00 UTC (rev 11908)
@@ -192,7 +192,7 @@
 	E->segment.Vx = (float *)malloc((E->lmesh.noz + 1) * sizeof(float));
 
 	E->F = (double *)malloc((E->mesh.nsd * E->lmesh.nnov + 1) * sizeof(double));
-	E->U = (double *)malloc((E->mesh.nsd * E->lmesh.nnov + 1) * sizeof(double));
+	E->U = (double *)malloc((E->lmesh.neq + 2) * sizeof(double));
 	E->T = (float *)malloc((E->lmesh.nno + 1) * sizeof(float));
 	E->C = (float *)malloc((E->lmesh.nno + 1) * sizeof(float));
 	E->CE = (float *)malloc((E->lmesh.nel + 1) * sizeof(float));
@@ -361,7 +361,7 @@
 
 	for(i = 1; i <= E->lmesh.nno; i++)
 		E->T[i] = E->buoyancy[i] = 0.0;
-	for(i = 0; i < E->lmesh.neq; i++)
+	for(i = 0; i <= E->lmesh.neq + 1; i++)
 		E->U[i] = 0.0;
 
 	set_up_nonmg_aliases(E);

Modified: mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c	2008-05-02 23:43:25 UTC (rev 11907)
+++ mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c	2008-05-02 23:56:00 UTC (rev 11908)
@@ -130,7 +130,7 @@
 		s1 = (double *)malloc((npno + 1) * sizeof(double));
 		s2 = (double *)malloc((npno + 1) * sizeof(double));
 		Ah = (double *)malloc((neq + 1) * sizeof(double));
-		u1 = (double *)malloc((neq + 1) * sizeof(double));
+		u1 = (double *)malloc((neq + 2) * sizeof(double));
 		been_here = 1;
 		timea = CPU_time0();
 	}