[cig-commits] r12257 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Tue Jun 17 13:31:30 PDT 2008
Author: becker
Date: 2008-06-17 13:31:30 -0700 (Tue, 17 Jun 2008)
New Revision: 12257
Modified:
mc/3D/CitcomS/trunk/lib/Drive_solvers.c
mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Size_does_matter.c
mc/3D/CitcomS/trunk/lib/Solver_conj_grad.c
mc/3D/CitcomS/trunk/lib/Solver_multigrid.c
mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c
mc/3D/CitcomS/trunk/lib/Topo_gravity.c
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Added a few more comments and made a few loops more efficient.
Modified: mc/3D/CitcomS/trunk/lib/Drive_solvers.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2008-06-17 19:15:51 UTC (rev 12256)
+++ mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2008-06-17 20:31:30 UTC (rev 12257)
@@ -93,6 +93,7 @@
solve_constrained_flow_iterative(E);
if (E->viscosity.SDEPV || E->viscosity.PDEPV) {
+ /* outer iterations for velocity dependent viscosity */
for (m=1;m<=E->sphere.caps_per_proc;m++) {
delta_U[m] = (double *)malloc(neq*sizeof(double));
Modified: mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/General_matrix_functions.c 2008-06-17 19:15:51 UTC (rev 12256)
+++ mc/3D/CitcomS/trunk/lib/General_matrix_functions.c 2008-06-17 20:31:30 UTC (rev 12257)
@@ -88,13 +88,15 @@
initial_time=CPU_time0();
if (!(E->control.NMULTIGRID || E->control.EMULTIGRID)) {
+ /* conjugate gradient solution */
+
cycles = E->control.v_steps_low;
residual = conj_grad(E,d0,F,acc,&cycles,high_lev);
valid = (residual < acc)? 1:0;
- }
+ } else {
+
+ /* solve using multigrid */
- else {
-
counts =0;
if(E->parallel.me==0){ /* output */
snprintf(message,200,"resi = %.6e for iter %d acc %.6e",residual,counts,acc);
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2008-06-17 19:15:51 UTC (rev 12256)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2008-06-17 20:31:30 UTC (rev 12257)
@@ -698,7 +698,8 @@
E->buoyancy[j] = (double *) malloc((nno+1)*sizeof(double));
E->gstress[j] = (float *) malloc((6*nno+1)*sizeof(float));
- E->stress[j] = (float *) malloc((12*nsf+1)*sizeof(float));
+ // TWB do we need this anymore XXX
+ //E->stress[j] = (float *) malloc((12*nsf+1)*sizeof(float));
for(i=1;i<=E->mesh.nsd;i++)
E->sphere.cap[j].TB[i] = (float *) malloc((nno+1)*sizeof(float));
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2008-06-17 19:15:51 UTC (rev 12256)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2008-06-17 20:31:30 UTC (rev 12257)
@@ -251,7 +251,8 @@
fprintf(fp2,"%3d %7d\n",j,E->lmesh.nsf);
for(i=1;i<=E->lmesh.nsf;i++) {
s = i*E->lmesh.noz;
- fprintf(fp2,"%.4e %.4e %.4e %.4e\n",topo[i],E->slice.shflux[j][i],E->sphere.cap[j].V[1][s],E->sphere.cap[j].V[2][s]);
+ fprintf(fp2,"%.4e %.4e %.4e %.4e\n",
+ topo[i],E->slice.shflux[j][i],E->sphere.cap[j].V[1][s],E->sphere.cap[j].V[2][s]);
}
}
fclose(fp2);
@@ -267,7 +268,8 @@
fprintf(fp2,"%3d %7d\n",j,E->lmesh.nsf);
for(i=1;i<=E->lmesh.nsf;i++) {
s = (i-1)*E->lmesh.noz + 1;
- fprintf(fp2,"%.4e %.4e %.4e %.4e\n",E->slice.tpgb[j][i],E->slice.bhflux[j][i],E->sphere.cap[j].V[1][s],E->sphere.cap[j].V[2][s]);
+ fprintf(fp2,"%.4e %.4e %.4e %.4e\n",
+ E->slice.tpgb[j][i],E->slice.bhflux[j][i],E->sphere.cap[j].V[1][s],E->sphere.cap[j].V[2][s]);
}
}
fclose(fp2);
Modified: mc/3D/CitcomS/trunk/lib/Size_does_matter.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Size_does_matter.c 2008-06-17 19:15:51 UTC (rev 12256)
+++ mc/3D/CitcomS/trunk/lib/Size_does_matter.c 2008-06-17 20:31:30 UTC (rev 12257)
@@ -195,6 +195,11 @@
}
+/*
+
+gets r,theta,phi coordinates at the integration points
+
+ */
void get_rtf_at_vpts(struct All_variables *E, int m, int lev, int el,
double rtf[4][9])
{
Modified: mc/3D/CitcomS/trunk/lib/Solver_conj_grad.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Solver_conj_grad.c 2008-06-17 19:15:51 UTC (rev 12256)
+++ mc/3D/CitcomS/trunk/lib/Solver_conj_grad.c 2008-06-17 20:31:30 UTC (rev 12257)
@@ -34,7 +34,7 @@
void solve_constrained_flow_iterative();
void cg_allocate_vars();
- E->build_forcing_term = assemble_forces_iterative;
+ E->build_forcing_term = assemble_forces_iterative;
E->solve_stokes_problem = solve_constrained_flow_iterative;
E->solver_allocate_vars = cg_allocate_vars;
Modified: mc/3D/CitcomS/trunk/lib/Solver_multigrid.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Solver_multigrid.c 2008-06-17 19:15:51 UTC (rev 12256)
+++ mc/3D/CitcomS/trunk/lib/Solver_multigrid.c 2008-06-17 20:31:30 UTC (rev 12257)
@@ -36,10 +36,11 @@
void solve_constrained_flow_iterative();
void mg_allocate_vars();
+ E->build_forcing_term = assemble_forces_iterative;
+ E->solve_stokes_problem = solve_constrained_flow_iterative;
E->solver_allocate_vars = mg_allocate_vars;
- E->build_forcing_term = assemble_forces_iterative;
- E->solve_stokes_problem = solve_constrained_flow_iterative;
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c 2008-06-17 19:15:51 UTC (rev 12256)
+++ mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c 2008-06-17 20:31:30 UTC (rev 12257)
@@ -105,9 +105,9 @@
double modified_plgndr_a();
int m,node,ll,mm,i,j,p;
- double t,f;
+ double t,f,mmf;
+
-
for(m=1;m<=E->sphere.caps_per_proc;m++) {
E->sphere.tablesplm[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
E->sphere.tablescosf[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
@@ -126,8 +126,9 @@
f=E->sx[m][2][node];
t=E->sx[m][1][node];
for (mm=0;mm<=E->output.llmax;mm++) {
- E->sphere.tablescosf[m][j][mm] = cos( (double)(mm)*f );
- E->sphere.tablessinf[m][j][mm] = sin( (double)(mm)*f );
+ mmf = (double)(mm)*f;
+ E->sphere.tablescosf[m][j][mm] = cos( mmf );
+ E->sphere.tablessinf[m][j][mm] = sin( mmf );
}
for (ll=0;ll<=E->output.llmax;ll++)
Modified: mc/3D/CitcomS/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2008-06-17 19:15:51 UTC (rev 12256)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2008-06-17 20:31:30 UTC (rev 12257)
@@ -41,7 +41,7 @@
void allocate_STD_mem();
void compute_nodal_stress();
void free_STD_mem();
- void get_surf_stress();
+ //void get_surf_stress();
int node,snode,m;
float *SXX[NCS],*SYY[NCS],*SXY[NCS],*SXZ[NCS],*SZY[NCS],*SZZ[NCS];
@@ -51,8 +51,9 @@
allocate_STD_mem(E, SXX, SYY, SZZ, SXY, SXZ, SZY, divv, vorv);
compute_nodal_stress(E, SXX, SYY, SZZ, SXY, SXZ, SZY, divv, vorv);
- if (E->parallel.me_loc[3]==E->parallel.nprocz-1)
- get_surf_stress(E,SXX,SYY,SZZ,SXY,SXZ,SZY);
+ // not needed ? TWB XXX
+ //if (E->parallel.me_loc[3]==E->parallel.nprocz-1)
+ //get_surf_stress(E,SXX,SYY,SZZ,SXY,SXZ,SZY);
topo_scaling1 = topo_scaling2 = 1.0;
@@ -60,9 +61,10 @@
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(snode=1;snode<=E->lmesh.nsf;snode++) {
node = E->surf_node[m][snode];
- tpg[m][snode] = -2*SZZ[m][node] + SZZ[m][node-1];
- tpgb[m][snode] = 2*SZZ[m][node-E->lmesh.noz+1]-SZZ[m][node-E->lmesh.noz+2];
- tpg[m][snode] = tpg[m][snode]*topo_scaling1;
+ tpg[m][snode] = -2*SZZ[m][node] + SZZ[m][node-1];
+ tpgb[m][snode] = 2*SZZ[m][node-E->lmesh.noz+1]- SZZ[m][node-E->lmesh.noz+2];
+
+ tpg[m][snode] = tpg[m][snode] *topo_scaling1;
tpgb[m][snode] = tpgb[m][snode]*topo_scaling2;
divg[m][snode] = 2*divv[m][node]-divv[m][node-1];
@@ -142,37 +144,37 @@
}
-void get_surf_stress(E,SXX,SYY,SZZ,SXY,SXZ,SZY)
- struct All_variables *E;
- float **SXX,**SYY,**SZZ,**SXY,**SXZ,**SZY;
-{
- int m,i,node,stride;
+/* void get_surf_stress(E,SXX,SYY,SZZ,SXY,SXZ,SZY) */
+/* struct All_variables *E; */
+/* float **SXX,**SYY,**SZZ,**SXY,**SXZ,**SZY; */
+/* { */
+/* int m,i,node,stride; */
- stride = E->lmesh.nsf*6;
+/* stride = E->lmesh.nsf*6; */
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for (node=1;node<=E->lmesh.nno;node++)
- if ( (node%E->lmesh.noz)==0 ) {
- i = node/E->lmesh.noz;
- E->stress[m][(i-1)*6+1] = SXX[m][node];
- E->stress[m][(i-1)*6+2] = SYY[m][node];
- E->stress[m][(i-1)*6+3] = SZZ[m][node];
- E->stress[m][(i-1)*6+4] = SXY[m][node];
- E->stress[m][(i-1)*6+5] = SXZ[m][node];
- E->stress[m][(i-1)*6+6] = SZY[m][node];
- }
- else if ( ((node+1)%E->lmesh.noz)==0 ) {
- i = (node+1)/E->lmesh.noz;
- E->stress[m][stride+(i-1)*6+1] = SXX[m][node];
- E->stress[m][stride+(i-1)*6+2] = SYY[m][node];
- E->stress[m][stride+(i-1)*6+3] = SZZ[m][node];
- E->stress[m][stride+(i-1)*6+4] = SXY[m][node];
- E->stress[m][stride+(i-1)*6+5] = SXZ[m][node];
- E->stress[m][stride+(i-1)*6+6] = SZY[m][node];
- }
+/* for(m=1;m<=E->sphere.caps_per_proc;m++) */
+/* for (node=1;node<=E->lmesh.nno;node++) */
+/* if ( (node%E->lmesh.noz)==0 ) { */
+/* i = node/E->lmesh.noz; */
+/* E->stress[m][(i-1)*6+1] = SXX[m][node]; */
+/* E->stress[m][(i-1)*6+2] = SYY[m][node]; */
+/* E->stress[m][(i-1)*6+3] = SZZ[m][node]; */
+/* E->stress[m][(i-1)*6+4] = SXY[m][node]; */
+/* E->stress[m][(i-1)*6+5] = SXZ[m][node]; */
+/* E->stress[m][(i-1)*6+6] = SZY[m][node]; */
+/* } */
+/* else if ( ((node+1)%E->lmesh.noz)==0 ) { */
+/* i = (node+1)/E->lmesh.noz; */
+/* E->stress[m][stride+(i-1)*6+1] = SXX[m][node]; */
+/* E->stress[m][stride+(i-1)*6+2] = SYY[m][node]; */
+/* E->stress[m][stride+(i-1)*6+3] = SZZ[m][node]; */
+/* E->stress[m][stride+(i-1)*6+4] = SXY[m][node]; */
+/* E->stress[m][stride+(i-1)*6+5] = SXZ[m][node]; */
+/* E->stress[m][stride+(i-1)*6+6] = SZY[m][node]; */
+/* } */
- return;
-}
+/* return; */
+/* } */
void compute_nodal_stress(struct All_variables *E,
@@ -189,7 +191,7 @@
float VV[4][9],Vxyz[9][9],Szz,Sxx,Syy,Sxy,Sxz,Szy,div,vor;
double dilation[9];
double pre[9],tww[9],rtf[4][9];
- double velo_scaling, stress_scaling;
+ double velo_scaling, stress_scaling, mass_fac;
struct Shape_function_dA *dOmega;
struct Shape_function_dx *GNx;
@@ -211,10 +213,14 @@
div = 0.0;
vor = 0.0;
- get_rtf_at_vpts(E, m, lev, e, rtf);
- velo_from_element(E,VV,m,e,sphere_key);
- dOmega = &(E->gDA[m][e]);
- GNx = &(E->gNX[m][e]);
+ get_rtf_at_vpts(E, m, lev, e, rtf);// gets r,theta,phi coordinates at the integration points
+ velo_from_element(E,VV,m,e,sphere_key); /* assign node-global
+ velocities to nodes
+ local to the
+ element */
+ dOmega = &(E->gDA[m][e]); /* Jacobian at integration points */
+ GNx = &(E->gNX[m][e]); /* derivatives of shape functions at
+ integration points */
/* Vxyz is the strain rate vector, whose relationship with
* the strain rate tensor (e) is that:
@@ -241,13 +247,17 @@
for(i=1;i<=ends;i++) {
tww[i] = 0.0;
- for(j=1;j<=vpts;j++)
+ for(j=1;j<=vpts;j++) /* weighting, consisting of Jacobian,
+ Gauss weight and shape function,
+ evaluated at integration points */
tww[i] += dOmega->vpt[j] * g_point[j].weight[E->mesh.nsd-1]
* E->N.vpt[GNVINDEX(i,j)];
}
- for(j=1;j<=vpts;j++) {
- for(i=1;i<=ends;i++) {
+ /* integrate over element */
+ for(j=1;j<=vpts;j++) { /* Gauss integration points */
+ for(i=1;i<=ends;i++) { /* nodes in element loop */
+ /* strain rate contributions from each node */
Vxyz[1][j]+=( VV[1][i]*GNx->vpt[GNVXINDEX(0,i,j)]
+ VV[3][i]*E->N.vpt[GNVINDEX(i,j)] )*rtf[3][j];
Vxyz[2][j]+=( (VV[2][i]*GNx->vpt[GNVXINDEX(1,i,j)]
@@ -273,22 +283,22 @@
}
}
- if(E->control.inv_gruneisen != 0) {
+ if(E->control.inv_gruneisen != 0) { /* isotropic component */
for(j=1;j<=vpts;j++)
dilation[j] = (Vxyz[1][j] + Vxyz[2][j] + Vxyz[3][j]) / 3.0;
}
for(j=1;j<=vpts;j++) {
- Sxx += 2.0 * pre[j] * (Vxyz[1][j] - dilation[j]);
+ Sxx += 2.0 * pre[j] * (Vxyz[1][j] - dilation[j]); /* */
Syy += 2.0 * pre[j] * (Vxyz[2][j] - dilation[j]);
Szz += 2.0 * pre[j] * (Vxyz[3][j] - dilation[j]);
- Sxy += pre[j] * Vxyz[4][j];
+ Sxy += pre[j] * Vxyz[4][j]; /* */
Sxz += pre[j] * Vxyz[5][j];
Szy += pre[j] * Vxyz[6][j];
- div += Vxyz[7][j]*dOmega->vpt[j];
- vor += Vxyz[8][j]*dOmega->vpt[j];
+ div += Vxyz[7][j]*dOmega->vpt[j]; /* divergence */
+ vor += Vxyz[8][j]*dOmega->vpt[j]; /* vorticity */
}
-
+ /* normalize by volume */
Sxx /= E->eco[m][e].area;
Syy /= E->eco[m][e].area;
Szz /= E->eco[m][e].area;
@@ -304,7 +314,7 @@
Syy -= E->P[m][e];
for(i=1;i<=ends;i++) {
- node = E->ien[m][e].node[i];
+ node = E->ien[m][e].node[i]; /* assign to global nodes */
SZZ[m][node] += tww[i] * Szz;
SXX[m][node] += tww[i] * Sxx;
SYY[m][node] += tww[i] * Syy;
@@ -331,14 +341,17 @@
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(node=1;node<=E->lmesh.nno;node++) {
- SZZ[m][node] = SZZ[m][node]*E->Mass[m][node]*stress_scaling;
- SXX[m][node] = SXX[m][node]*E->Mass[m][node]*stress_scaling;
- SYY[m][node] = SYY[m][node]*E->Mass[m][node]*stress_scaling;
- SXY[m][node] = SXY[m][node]*E->Mass[m][node]*stress_scaling;
- SXZ[m][node] = SXZ[m][node]*E->Mass[m][node]*stress_scaling;
- SZY[m][node] = SZY[m][node]*E->Mass[m][node]*stress_scaling;
- vorv[m][node] = vorv[m][node]*E->Mass[m][node]*velo_scaling;
- divv[m][node] = divv[m][node]*E->Mass[m][node]*velo_scaling;
+ mass_fac = E->Mass[m][node]*stress_scaling;
+ SZZ[m][node] *= mass_fac;
+ SXX[m][node] *= mass_fac;
+ SYY[m][node] *= mass_fac;
+ SXY[m][node] *= mass_fac;
+ SXZ[m][node] *= mass_fac;
+ SZY[m][node] *= mass_fac;
+
+ mass_fac = E->Mass[m][node]*velo_scaling;
+ vorv[m][node] *= mass_fac;
+ divv[m][node] *= mass_fac;
}
/* assign stress to all the nodes */
@@ -375,8 +388,11 @@
for(i=1; i<=E->lmesh.noy; i++)
for(j=1; j<=E->lmesh.nox; j++)
for(k=1; k<=E->lmesh.noz; k++) {
+
n = k+(j-1)*E->lmesh.noz+(i-1)*E->lmesh.nox*E->lmesh.noz;
+
for(d=1; d<=E->mesh.nsd; d++)
+
if(E->node[m][n] & sbc_flag[d]) {
if(i==1)
E->gstress[m][(n-1)*6+stress_index[d][2]] = E->sbc.SB[m][SIDE_WEST][d][ E->sbc.node[m][n] ];
@@ -428,12 +444,16 @@
* and dimensionalized (data.density). dlayer needs to be dimensionalized.
*/
- int m,k,ll,mm,node,i,j,p,noz,snode;
- float *TT[NCS],radius,*geoid[2],dlayer,con1,grav,scaling2,scaling;
+ int m,k,ll,mm,node,i,j,p,noz,snode,nxnz;
+ float *TT[NCS],radius,*geoid[2],dlayer,con1,grav,scaling2,scaling,radius_m;
double buoy2rho;
void sphere_expansion();
void sum_across_depth_sph1();
+ /* some constants */
+ nxnz = E->lmesh.nox*E->lmesh.noz;
+ radius_m = E->data.radius_km*1e3;
+
/* scale for buoyancy */
scaling2 = -E->data.therm_exp*E->data.ref_temperature*E->data.density
/ E->control.Atemp;
@@ -464,7 +484,7 @@
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.noy;i++)
for(j=1;j<=E->lmesh.nox;j++) {
- node=k+(j-1)*E->lmesh.noz+(i-1)*E->lmesh.nox*E->lmesh.noz;
+ node= k + (j-1)*E->lmesh.noz + (i-1)*nxnz;
p = j + (i-1)*E->lmesh.nox;
TT[m][p] = (E->buoyancy[m][node]+E->buoyancy[m][node+1])
* 0.5 * buoy2rho;
@@ -474,15 +494,14 @@
sphere_expansion(E,TT,geoid[0],geoid[1]);
/* thickness of the layer */
- dlayer = (E->sx[1][3][k+1]-E->sx[1][3][k])*E->data.radius_km*1e3;
+ dlayer = (E->sx[1][3][k+1]-E->sx[1][3][k])*radius_m;
/* mean radius of the layer */
radius = (E->sx[1][3][k+1]+E->sx[1][3][k])*0.5;
/* geoid contribution of density at this layer, ignore degree-0 term */
for (ll=1;ll<=E->output.llmax;ll++) {
- con1 = scaling * dlayer * pow(radius,((double)(ll+2)))
- / (2.0*ll+1.0);
+ con1 = scaling * dlayer * pow(radius,((double)(ll+2))) / (2.0*ll+1.0);
for (mm=0;mm<=ll;mm++) {
p = E->sphere.hindex[ll][mm];
harm_geoid[0][p] += con1*geoid[0][p];
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2008-06-17 19:15:51 UTC (rev 12256)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2008-06-17 20:31:30 UTC (rev 12257)
@@ -752,7 +752,7 @@
double *surf_det[NCS][5];
double *SinCos[MAX_LEVELS][NCS][4];
- float *stress[NCS];
+ //float *stress[NCS];
float *gstress[NCS];
float *Fas670[NCS],*Fas410[NCS],*Fas670_b[NCS],*Fas410_b[NCS];
float *Fascmb[NCS],*Fascmb_b[NCS];
More information about the cig-commits
mailing list