[cig-commits] r7936 - mc/3D/CitcomS/trunk/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Sep 5 15:35:31 PDT 2007
Author: tan2
Date: 2007-09-05 15:35:30 -0700 (Wed, 05 Sep 2007)
New Revision: 7936
Modified:
mc/3D/CitcomS/trunk/lib/Checkpoints.c
mc/3D/CitcomS/trunk/lib/Construct_arrays.c
mc/3D/CitcomS/trunk/lib/Drive_solvers.c
mc/3D/CitcomS/trunk/lib/Element_calculations.c
mc/3D/CitcomS/trunk/lib/Full_obsolete.c
mc/3D/CitcomS/trunk/lib/Full_parallel_related.c
mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Obsolete.c
mc/3D/CitcomS/trunk/lib/Regional_obsolete.c
mc/3D/CitcomS/trunk/lib/Regional_parallel_related.c
mc/3D/CitcomS/trunk/lib/Solver_multigrid.c
mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
Log:
Fixed a few memory problems.
Shrink the size of several arrays related to Stokes eqn.
Moved jacobi() to Obsolete.c
Modified: mc/3D/CitcomS/trunk/lib/Checkpoints.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Checkpoints.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/Checkpoints.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -454,7 +454,7 @@
fwrite(E->NP[m], sizeof(float), E->lmesh.nno+1, fp);
/* velocity at equation points */
- fwrite(E->U[m], sizeof(double), E->lmesh.neq+2, fp);
+ fwrite(E->U[m], sizeof(double), E->lmesh.neq, fp);
/* viscosity at quadrature points and node points */
fwrite(E->EVI[lev][m], sizeof(float),
@@ -484,7 +484,7 @@
fread(E->NP[m], sizeof(float), E->lmesh.nno+1, fp);
/* velocity at equation points */
- fread(E->U[m], sizeof(double), E->lmesh.neq+2, fp);
+ fread(E->U[m], sizeof(double), E->lmesh.neq, fp);
/* viscosity at quadrature points and node points */
fread(E->EVI[lev][m], sizeof(float),
Modified: mc/3D/CitcomS/trunk/lib/Construct_arrays.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Construct_arrays.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/Construct_arrays.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -254,12 +254,12 @@
noy = E->lmesh.NOY[lev];
nox = E->lmesh.NOX[lev];
max_eqn = 14*dims;
- matrix = max_eqn*(nno+3);
+ matrix = max_eqn*nno;
- E->Node_map[lev][m]=(int *) malloc ((matrix+3)*sizeof(int));
+ E->Node_map[lev][m]=(int *) malloc (matrix*sizeof(int));
- for(i=0;i<=matrix;i++)
- E->Node_map[lev][m][i] = neq+1; /* DANGER !!! */
+ for(i=0;i<matrix;i++)
+ E->Node_map[lev][m][i] = neq; /* neq indicates an invalid eqn # */
for (ii=1;ii<=noy;ii++)
for (jj=1;jj<=nox;jj++)
@@ -290,13 +290,21 @@
}
}
- E->Eqn_k1[lev][m] = (higher_precision *)malloc((matrix+5)*sizeof(higher_precision));
- E->Eqn_k2[lev][m] = (higher_precision *)malloc((matrix+5)*sizeof(higher_precision));
- E->Eqn_k3[lev][m] = (higher_precision *)malloc((matrix+5)*sizeof(higher_precision));
+ E->Eqn_k1[lev][m] = (higher_precision *)malloc(matrix*sizeof(higher_precision));
+ E->Eqn_k2[lev][m] = (higher_precision *)malloc(matrix*sizeof(higher_precision));
+ E->Eqn_k3[lev][m] = (higher_precision *)malloc(matrix*sizeof(higher_precision));
- E->mesh.matrix_size[lev] = matrix + 1;
- } /* end for level and m */
+ E->mesh.matrix_size[lev] = matrix;
+ if(E->control.verbose) {
+ fprintf(E->fp_out, "output Node_map lev=%d m=%d\n", lev, m);
+ fprintf(E->fp_out, "neq=%d nno=%d max_eqn=%d matrix=%d\n", neq, nno, max_eqn, matrix);
+ for(i=0;i<matrix;i++)
+ fprintf(E->fp_out, "%d %d\n", i, E->Node_map[lev][m][i]);
+ }
+
+ } /* end for level and m */
+
return;
}
@@ -332,9 +340,9 @@
neq=E->lmesh.NEQ[level];
nel=E->lmesh.NEL[level];
nno=E->lmesh.NNO[level];
- for(i=0;i<=(neq+1);i++)
+ for(i=0;i<neq;i++)
E->BI[level][m][i] = zero;
- for(i=0;i<=E->mesh.matrix_size[level];i++) {
+ for(i=0;i<E->mesh.matrix_size[level];i++) {
E->Eqn_k1[level][m][i] = zero;
E->Eqn_k2[level][m][i] = zero;
E->Eqn_k3[level][m][i] = zero;
@@ -464,7 +472,7 @@
for(level=E->mesh.gridmax;level>=E->mesh.gridmin;level--) {
for (m=1;m<=E->sphere.caps_per_proc;m++) {
- for(j=0;j<E->lmesh.NEQ[level];j++)
+ for(j=0;j<=E->lmesh.NEQ[level];j++)
E->temp[m][j]=0.0;
for(i=1;i<=E->lmesh.NNO[level];i++) {
Modified: mc/3D/CitcomS/trunk/lib/Drive_solvers.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -94,9 +94,9 @@
if (E->viscosity.SDEPV || E->viscosity.PDEPV) {
for (m=1;m<=E->sphere.caps_per_proc;m++) {
- delta_U[m] = (double *)malloc((neq+2)*sizeof(double));
- oldU[m] = (double *)malloc((neq+2)*sizeof(double));
- for(i=0;i<=neq;i++)
+ delta_U[m] = (double *)malloc(neq*sizeof(double));
+ oldU[m] = (double *)malloc(neq*sizeof(double));
+ for(i=0;i<neq;i++)
oldU[m][i]=0.0;
}
@@ -185,9 +185,9 @@
if (E->viscosity.SDEPV || E->viscosity.PDEPV) {
for (m=1;m<=E->sphere.caps_per_proc;m++) {
- delta_U[m] = (double *)malloc((neq+2)*sizeof(double));
- oldU[m] = (double *)malloc((neq+2)*sizeof(double));
- for(i=0;i<=neq;i++)
+ delta_U[m] = (double *)malloc(neq*sizeof(double));
+ oldU[m] = (double *)malloc(neq*sizeof(double));
+ for(i=0;i<neq;i++)
oldU[m][i]=0.0;
}
Modified: mc/3D/CitcomS/trunk/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Element_calculations.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/Element_calculations.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -525,7 +525,7 @@
for(e=0;e<=neq;e++)
Au[m][e]=0.0;
- u[m][neq]=0.0;
+ u[m][neq] = 0.0;
for(e=1;e<=nno;e++) {
@@ -543,15 +543,12 @@
B3=E->Eqn_k3[level][m]+(e-1)*max_eqn;
for(i=3;i<max_eqn;i++) {
- if (C[i] != neq+1) {
UU = u[m][C[i]];
Au[m][eqn1] += B1[i]*UU;
Au[m][eqn2] += B2[i]*UU;
Au[m][eqn3] += B3[i]*UU;
- }
}
for(i=0;i<max_eqn;i++)
- if (C[i] != neq+1)
Au[m][C[i]] += B1[i]*U1+B2[i]*U2+B3[i]*U3;
} /* end for e */
@@ -740,7 +737,7 @@
nel=E->lmesh.NEL[lev];
neq=E->lmesh.NEQ[lev];
- for(i=0;i<=neq;i++)
+ for(i=0;i<neq;i++)
gradP[m][i] = 0.0;
for(e=1;e<=nel;e++) {
Modified: mc/3D/CitcomS/trunk/lib/Full_obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_obsolete.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/Full_obsolete.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -93,7 +93,7 @@
processors = (int *)malloc((E->parallel.nprocz+2)*sizeof(int));
- SD = (double *)malloc((E->lmesh.NEQ[lev]+2)*sizeof(double));
+ SD = (double *)malloc((E->lmesh.NEQ[lev])*sizeof(double));
rootid = E->parallel.me_sph*E->parallel.nprocz; /* which is the bottom cpu */
@@ -127,7 +127,7 @@
MPI_Send(SD,E->lmesh.NEQ[lev],MPI_DOUBLE,processors[d],rootid,E->parallel.world);
}
else
- for (i=0;i<=E->lmesh.NEQ[lev];i++)
+ for (i=0;i<E->lmesh.NEQ[lev];i++)
AUo[m][i] = SD[i];
}
else
@@ -164,7 +164,7 @@
processors = (int *)malloc((E->parallel.nprocz+2)*sizeof(int));
- RV = (double *)malloc((E->lmesh.NEQ[lev]+2)*sizeof(double));
+ RV = (double *)malloc((E->lmesh.NEQ[lev])*sizeof(double));
rootid = E->parallel.me_sph*E->parallel.nprocz; /* which is the bottom cpu */
Modified: mc/3D/CitcomS/trunk/lib/Full_parallel_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_parallel_related.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/Full_parallel_related.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -486,7 +486,7 @@
fprintf(E->fp_out,"ii=%d %d %d \n",ii,i,E->parallel.NODE[lev][m][i].bound[ii]);
lnode=0;
- for (node=1;node<=E->lmesh.nno;node++)
+ for (node=1;node<=E->lmesh.NNO[lev];node++)
if((E->NODE[lev][m][node] & SKIP)) {
lnode++;
fprintf(E->fp_out,"skip %d %d \n",lnode,node);
Modified: mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/General_matrix_functions.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/General_matrix_functions.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -199,7 +199,6 @@
void n_assemble_del2_u();
void strip_bcs_from_residual();
void gauss_seidel();
- void jacobi();
double conj_grad();
double multi_grid();
@@ -305,7 +304,6 @@
void project_vector();
int m,i,j,Vn,Vnmax,cycles;
double alpha,beta;
- void jacobi();
void gauss_seidel();
void element_gauss_seidel();
void e_assemble_del2_u();
@@ -330,12 +328,12 @@
for(i=E->mesh.levmin;i<=E->mesh.levmax;i++)
for(m=1;m<=E->sphere.caps_per_proc;m++) {
- del_vel[i][m]=(double *)malloc((E->lmesh.NEQ[i] + 2)*sizeof(double));
- AU[i][m] = (double *)malloc((E->lmesh.NEQ[i] + 2)*sizeof(double));
- vel[i][m]=(double *)malloc((E->lmesh.NEQ[i]+2)*sizeof(double));
- res[i][m]=(double *)malloc((E->lmesh.NEQ[i]+2)*sizeof(double));
+ del_vel[i][m]=(double *)malloc((E->lmesh.NEQ[i]+1)*sizeof(double));
+ AU[i][m] = (double *)malloc((E->lmesh.NEQ[i]+1)*sizeof(double));
+ vel[i][m]=(double *)malloc((E->lmesh.NEQ[i]+1)*sizeof(double));
+ res[i][m]=(double *)malloc((E->lmesh.NEQ[i])*sizeof(double));
if (i<E->mesh.levmax)
- fl[i][m]=(double *)malloc((E->lmesh.NEQ[i]+2)*sizeof(double));
+ fl[i][m]=(double *)malloc((E->lmesh.NEQ[i])*sizeof(double));
}
Vnmax = E->control.mg_cycle;
@@ -479,11 +477,11 @@
steps = *cycles;
for(m=1;m<=E->sphere.caps_per_proc;m++) {
- r0[m] = (double *)malloc((1+E->lmesh.NEQ[mem_lev])*sizeof(double));
- r1[m] = (double *)malloc((1+E->lmesh.NEQ[mem_lev])*sizeof(double));
- r2[m] = (double *)malloc((1+E->lmesh.NEQ[mem_lev])*sizeof(double));
- z0[m] = (double *)malloc((1+E->lmesh.NEQ[mem_lev])*sizeof(double));
- z1[m] = (double *)malloc((1+E->lmesh.NEQ[mem_lev])*sizeof(double));
+ r0[m] = (double *)malloc(E->lmesh.NEQ[mem_lev]*sizeof(double));
+ r1[m] = (double *)malloc(E->lmesh.NEQ[mem_lev]*sizeof(double));
+ r2[m] = (double *)malloc(E->lmesh.NEQ[mem_lev]*sizeof(double));
+ z0[m] = (double *)malloc(E->lmesh.NEQ[mem_lev]*sizeof(double));
+ z1[m] = (double *)malloc(E->lmesh.NEQ[mem_lev]*sizeof(double));
p1[m] = (double *)malloc((1+E->lmesh.NEQ[mem_lev])*sizeof(double));
p2[m] = (double *)malloc((1+E->lmesh.NEQ[mem_lev])*sizeof(double));
Ap[m] = (double *)malloc((1+E->lmesh.NEQ[mem_lev])*sizeof(double));
@@ -611,7 +609,7 @@
steps=*cycles;
for (m=1;m<=E->sphere.caps_per_proc;m++) {
- dd[m] = (double *)malloc((neq+1)*sizeof(double));
+ dd[m] = (double *)malloc(neq*sizeof(double));
vis[m] = (int *)malloc((nno+1)*sizeof(int));
}
elt_k=(double *)malloc((24*24)*sizeof(double));
@@ -736,126 +734,7 @@
return;
}
-void jacobi(E,d0,F,Ad,acc,cycles,level,guess)
- struct All_variables *E;
- double **d0;
- double **F,**Ad;
- double acc;
- int *cycles;
- int level;
- int guess;
-{
- int count,i,j,k,l,m,ns,steps;
- int *C;
- int eqn1,eqn2,eqn3,gneq;
-
- void n_assemble_del2_u();
-
- double sum1,sum2,sum3,residual,global_vdot(),U1,U2,U3;
-
- double *r1[NCS];
-
- higher_precision *B1,*B2,*B3;
-
-
- const int dims=E->mesh.nsd;
- const int ends=enodes[dims];
- const int n=loc_mat_size[E->mesh.nsd];
- const int neq=E->lmesh.NEQ[level];
- const int num_nodes=E->lmesh.NNO[level];
- const int nox=E->lmesh.NOX[level];
- const int noz=E->lmesh.NOY[level];
- const int noy=E->lmesh.NOZ[level];
- const int max_eqn=14*dims;
-
- gneq = E->mesh.NEQ[level];
-
- steps=*cycles;
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- r1[m] = (double *)malloc((E->lmesh.neq+1)*sizeof(double));
-
- if(guess) {
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- d0[m][neq]=0.0;
- n_assemble_del2_u(E,d0,Ad,level,1);
- }
- else
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<neq;i++) {
- d0[m][i]=Ad[m][i]=0.0;
- }
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<neq;i++)
- r1[m][i]=F[m][i]-Ad[m][i];
-
-
- count = 0;
-
- while (count < steps) {
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.NNO[level];i++) {
- eqn1=E->ID[level][m][i].doff[1];
- eqn2=E->ID[level][m][i].doff[2];
- eqn3=E->ID[level][m][i].doff[3];
- d0[m][eqn1] += r1[m][eqn1]*E->BI[level][m][eqn1];
- d0[m][eqn2] += r1[m][eqn2]*E->BI[level][m][eqn2];
- d0[m][eqn3] += r1[m][eqn3]*E->BI[level][m][eqn3];
- }
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<neq;i++)
- Ad[m][i]=0.0;
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.NNO[level];i++) {
- eqn1=E->ID[level][m][i].doff[1];
- eqn2=E->ID[level][m][i].doff[2];
- eqn3=E->ID[level][m][i].doff[3];
- U1 = d0[m][eqn1];
- U2 = d0[m][eqn2];
- U3 = d0[m][eqn3];
-
- C=E->Node_map[level][m]+(i-1)*max_eqn;
- B1=E->Eqn_k1[level][m]+(i-1)*max_eqn;
- B2=E->Eqn_k2[level][m]+(i-1)*max_eqn;
- B3=E->Eqn_k3[level][m]+(i-1)*max_eqn;
-
- for(j=3;j<max_eqn;j++) {
- Ad[m][eqn1] += B1[j]*d0[m][C[j]];
- Ad[m][eqn2] += B2[j]*d0[m][C[j]];
- Ad[m][eqn3] += B3[j]*d0[m][C[j]];
- }
-
- for(j=0;j<max_eqn;j++) {
- Ad[m][C[j]] += B1[j]*U1 + B2[j]*U2 + B3[j]*U3;
- }
- } /* end for i and m */
-
- (E->solver.exchange_id_d)(E, Ad, level);
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<neq;i++)
- r1[m][i] = F[m][i] - Ad[m][i];
-
- /* residual = sqrt(global_vdot(E,r1,r1,level))/gneq;
-
- if(E->parallel.me==0)fprintf(stderr,"residuall =%.5e for %d\n",residual,count);
-*/ count++;
- }
-
- *cycles=count;
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- free((double*) r1[m]);
-
- return;
-
- }
-
-
/* ============================================================================
Multigrid Gauss-Seidel relaxation scheme which requires the storage of local
information, otherwise some other method is required. NOTE this is a bit worse
@@ -903,13 +782,11 @@
sor = 1.3;
if(guess) {
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- d0[m][neq]=0.0;
n_assemble_del2_u(E,d0,Ad,level,1);
}
else
for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<=neq;i++) {
+ for(i=0;i<neq;i++) {
d0[m][i]=Ad[m][i]=zeroo;
}
@@ -922,6 +799,9 @@
E->temp[m][j] = zeroo;
for (m=1;m<=E->sphere.caps_per_proc;m++)
+ Ad[m][neq] = zeroo;
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.NNO[level];i++)
if(E->NODE[level][m][i] & OFFSIDE) {
@@ -1010,8 +890,6 @@
*/
}
- /* parallel_process_termination();
-*/
*cycles=count;
return;
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -725,10 +725,10 @@
E->lmesh.neq = E->lmesh.nnov * E->mesh.nsd;
E->temp[j] = (double *) malloc((E->lmesh.neq+1)*sizeof(double));
- E->temp1[j] = (double *) malloc((E->lmesh.neq+1)*sizeof(double));
- E->F[j] = (double *) malloc((E->lmesh.neq+1)*sizeof(double));
- E->U[j] = (double *) malloc((E->lmesh.neq+2)*sizeof(double));
- E->u1[j] = (double *) malloc((E->lmesh.neq+2)*sizeof(double));
+ E->temp1[j] = (double *) malloc(E->lmesh.neq*sizeof(double));
+ E->F[j] = (double *) malloc(E->lmesh.neq*sizeof(double));
+ E->U[j] = (double *) malloc((E->lmesh.neq+1)*sizeof(double));
+ E->u1[j] = (double *) malloc((E->lmesh.neq+1)*sizeof(double));
for(i=1;i<=E->mesh.nsd;i++) {
@@ -737,7 +737,7 @@
E->sphere.cap[j].Vprev[i] = (float *) malloc((E->lmesh.nnov+1)*sizeof(float));
}
- for(i=0;i<=E->lmesh.neq;i++)
+ for(i=0;i<E->lmesh.neq;i++)
E->U[j][i] = E->temp[j][i] = E->temp1[j][i] = 0.0;
if(E->control.tracer == 1) {
@@ -763,13 +763,13 @@
for (j=1;j<=E->sphere.caps_per_proc;j++) {
E->lmesh.NEQ[l] = E->lmesh.NNOV[l] * E->mesh.nsd;
- E->BI[l][j] = (double *) malloc((E->lmesh.NEQ[l]+2)*sizeof(double));
+ E->BI[l][j] = (double *) malloc((E->lmesh.NEQ[l])*sizeof(double));
k = (E->lmesh.NOX[l]*E->lmesh.NOZ[l]+E->lmesh.NOX[l]*E->lmesh.NOY[l]+
E->lmesh.NOY[l]*E->lmesh.NOZ[l])*6;
E->zero_resid[l][j] = (int *) malloc((k+2)*sizeof(int));
E->parallel.Skip_id[l][j] = (int *) malloc((k+2)*sizeof(int));
- for(i=0;i<E->lmesh.NEQ[l]+2;i++) {
+ for(i=0;i<E->lmesh.NEQ[l];i++) {
E->BI[l][j][i]=0.0;
}
@@ -1153,6 +1153,9 @@
/* break if no more field */
if(prev == NULL) break;
+ /* skip if empty */
+ if(prev[0] == '\0') continue;
+
/* strip off leading and trailing whitespaces */
prev = strip(prev);
Modified: mc/3D/CitcomS/trunk/lib/Obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Obsolete.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/Obsolete.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -965,6 +965,130 @@
}
+/* ========================================================== */
+/* From General_matrix_functions.c */
+/* =========================================================== */
+
+void jacobi(E,d0,F,Ad,acc,cycles,level,guess)
+ struct All_variables *E;
+ double **d0;
+ double **F,**Ad;
+ double acc;
+ int *cycles;
+ int level;
+ int guess;
+{
+
+ int count,i,j,k,l,m,ns,steps;
+ int *C;
+ int eqn1,eqn2,eqn3,gneq;
+
+ void n_assemble_del2_u();
+
+ double sum1,sum2,sum3,residual,global_vdot(),U1,U2,U3;
+
+ double *r1[NCS];
+
+ higher_precision *B1,*B2,*B3;
+
+
+ const int dims=E->mesh.nsd;
+ const int ends=enodes[dims];
+ const int n=loc_mat_size[E->mesh.nsd];
+ const int neq=E->lmesh.NEQ[level];
+ const int num_nodes=E->lmesh.NNO[level];
+ const int nox=E->lmesh.NOX[level];
+ const int noz=E->lmesh.NOY[level];
+ const int noy=E->lmesh.NOZ[level];
+ const int max_eqn=14*dims;
+
+ gneq = E->mesh.NEQ[level];
+
+ steps=*cycles;
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ r1[m] = (double *)malloc(E->lmesh.neq*sizeof(double));
+
+ if(guess) {
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ d0[m][neq]=0.0;
+ n_assemble_del2_u(E,d0,Ad,level,1);
+ }
+ else
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=0;i<=neq;i++) {
+ d0[m][i]=Ad[m][i]=0.0;
+ }
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=0;i<neq;i++)
+ r1[m][i]=F[m][i]-Ad[m][i];
+
+
+ count = 0;
+
+ while (count < steps) {
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=E->lmesh.NNO[level];i++) {
+ eqn1=E->ID[level][m][i].doff[1];
+ eqn2=E->ID[level][m][i].doff[2];
+ eqn3=E->ID[level][m][i].doff[3];
+ d0[m][eqn1] += r1[m][eqn1]*E->BI[level][m][eqn1];
+ d0[m][eqn2] += r1[m][eqn2]*E->BI[level][m][eqn2];
+ d0[m][eqn3] += r1[m][eqn3]*E->BI[level][m][eqn3];
+ }
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=0;i<=neq;i++)
+ Ad[m][i]=0.0;
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=E->lmesh.NNO[level];i++) {
+ eqn1=E->ID[level][m][i].doff[1];
+ eqn2=E->ID[level][m][i].doff[2];
+ eqn3=E->ID[level][m][i].doff[3];
+ U1 = d0[m][eqn1];
+ U2 = d0[m][eqn2];
+ U3 = d0[m][eqn3];
+
+ C=E->Node_map[level][m]+(i-1)*max_eqn;
+ B1=E->Eqn_k1[level][m]+(i-1)*max_eqn;
+ B2=E->Eqn_k2[level][m]+(i-1)*max_eqn;
+ B3=E->Eqn_k3[level][m]+(i-1)*max_eqn;
+
+ for(j=3;j<max_eqn;j++) {
+ Ad[m][eqn1] += B1[j]*d0[m][C[j]];
+ Ad[m][eqn2] += B2[j]*d0[m][C[j]];
+ Ad[m][eqn3] += B3[j]*d0[m][C[j]];
+ }
+
+ for(j=0;j<max_eqn;j++) {
+ Ad[m][C[j]] += B1[j]*U1 + B2[j]*U2 + B3[j]*U3;
+ }
+ } /* end for i and m */
+
+ (E->solver.exchange_id_d)(E, Ad, level);
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=0;i<neq;i++)
+ r1[m][i] = F[m][i] - Ad[m][i];
+
+ /* residual = sqrt(global_vdot(E,r1,r1,level))/gneq;
+
+ if(E->parallel.me==0)fprintf(stderr,"residuall =%.5e for %d\n",residual,count);
+*/ count++;
+ }
+
+ *cycles=count;
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ free((double*) r1[m]);
+
+ return;
+
+ }
+
+
/* version */
/* $Id$ */
Modified: mc/3D/CitcomS/trunk/lib/Regional_obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_obsolete.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/Regional_obsolete.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -96,7 +96,7 @@
processors = (int *)malloc((E->parallel.nprocz+2)*sizeof(int));
- SD = (double *)malloc((E->lmesh.NEQ[lev]+2)*sizeof(double));
+ SD = (double *)malloc((E->lmesh.NEQ[lev])*sizeof(double));
rootid = E->parallel.me_sph*E->parallel.nprocz; /* which is the bottom cpu */
@@ -168,7 +168,7 @@
processors = (int *)malloc((E->parallel.nprocz+2)*sizeof(int));
- RV = (double *)malloc((E->lmesh.NEQ[lev]+2)*sizeof(double));
+ RV = (double *)malloc((E->lmesh.NEQ[lev])*sizeof(double));
rootid = E->parallel.me_sph*E->parallel.nprocz; /* which is the bottom cpu */
Modified: mc/3D/CitcomS/trunk/lib/Regional_parallel_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_parallel_related.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/Regional_parallel_related.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -434,7 +434,7 @@
fprintf(E->fp_out,"ii=%d %d %d \n",ii,i,E->parallel.NODE[lev][m][i].bound[ii]);
lnode=0;
- for (node=1;node<=E->lmesh.nno;node++)
+ for (node=1;node<=E->lmesh.NNO[lev];node++)
if((E->NODE[lev][m][node] & SKIP)) {
lnode++;
fprintf(E->fp_out,"skip %d %d \n",lnode,node);
Modified: mc/3D/CitcomS/trunk/lib/Solver_multigrid.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Solver_multigrid.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/Solver_multigrid.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -149,7 +149,7 @@
assert(start_lev != E->mesh.levmax /* un_injection */);
for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=neq;i++)
+ for(i=1;i<neq;i++)
AU[m][i]=0.0;
for(m=1;m<=E->sphere.caps_per_proc;m++)
@@ -460,7 +460,7 @@
from_rtf_to_xyz(E,start_lev,AU,E->temp);
for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<=neq_minus;i++)
+ for(i=0;i<neq_minus;i++)
E->temp1[m][i] = 0.0;
/* smooth in xyz coordinates */
Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-09-05 22:26:43 UTC (rev 7935)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-09-05 22:35:30 UTC (rev 7936)
@@ -160,7 +160,7 @@
lev = E->mesh.levmax;
for (m=1; m<=E->sphere.caps_per_proc; m++) {
- F[m] = (double *)malloc((neq+1)*sizeof(double));
+ F[m] = (double *)malloc(neq*sizeof(double));
r1[m] = (double *)malloc((npno+1)*sizeof(double));
r2[m] = (double *)malloc((npno+1)*sizeof(double));
z1[m] = (double *)malloc((npno+1)*sizeof(double));
@@ -398,7 +398,7 @@
st[m] = (double *)malloc((npno+1)*sizeof(double));
t0[m] = (double *)malloc((npno+1)*sizeof(double));
- u0[m] = (double *)malloc((neq+1)*sizeof(double));
+ u0[m] = (double *)malloc(neq*sizeof(double));
}
time0 = CPU_time0();
@@ -629,8 +629,8 @@
double global_vdot(),global_pdot();
for (m=1;m<=E->sphere.caps_per_proc;m++) {
- old_v[m] = (double *)malloc((neq+1)*sizeof(double));
- diff_v[m] = (double *)malloc((neq+1)*sizeof(double));
+ old_v[m] = (double *)malloc(neq*sizeof(double));
+ diff_v[m] = (double *)malloc(neq*sizeof(double));
old_p[m] = (double *)malloc((npno+1)*sizeof(double));
diff_p[m] = (double *)malloc((npno+1)*sizeof(double));
}
More information about the cig-commits
mailing list