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

[tan2 at geodynamics.org]

Author: tan2
Date: 2007-09-12 12:24:15 -0700 (Wed, 12 Sep 2007)
New Revision: 7956

Modified:
   mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
Log:
Keep the load vector intact in the bicg solver

Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c	2007-09-12 19:23:56 UTC (rev 7955)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c	2007-09-12 19:24:15 UTC (rev 7956)
@@ -351,7 +351,7 @@
  */
 
 static float solve_Ahat_p_fhat_BiCG(struct All_variables *E,
-                                    double **V, double **P, double **F,
+                                    double **V, double **P, double **FF,
                                     double imp, int *steps_max)
 {
     void assemble_div_rho_u();
@@ -373,6 +373,7 @@
     double r0dotrt, r1dotrt;
     double residual, dpressure, dvelocity;
 
+    double *F[NCS];
     double *r1[NCS], *r2[NCS], *pt[NCS], *p1[NCS], *p2[NCS];
     double *rt[NCS], *v0[NCS], *s0[NCS], *st[NCS], *t0[NCS];
     double *u0[NCS];
@@ -387,6 +388,7 @@
     lev = E->mesh.levmax;
 
     for (m=1; m<=E->sphere.caps_per_proc; m++)   {
+        F[m] = (double *)malloc(neq*sizeof(double));
         r1[m] = (double *)malloc((npno+1)*sizeof(double));
         r2[m] = (double *)malloc((npno+1)*sizeof(double));
         pt[m] = (double *)malloc((npno+1)*sizeof(double));
@@ -404,6 +406,13 @@
     time0 = CPU_time0();
     count = 0;
 
+    /* copy the original force vector since we need to keep it intact
+       between iterations */
+    for(m=1;m<=E->sphere.caps_per_proc;m++)
+        for(j=0;j<neq;j++)
+            F[m][j] = FF[m][j];
+
+
     /* calculate the initial velocity residual */
     v_res = initial_vel_residual(E, V, P, F, imp);
 
@@ -587,6 +596,7 @@
 
 
     for(m=1; m<=E->sphere.caps_per_proc; m++) {
+    	free((void *) F[m]);
         free((void *) r1[m]);
         free((void *) r2[m]);
         free((void *) pt[m]);