[cig-commits] r18376 - in mc/1D/hc/trunk: . anisotropic_viscosity
becker at geodynamics.org
becker at geodynamics.org
Sun May 15 15:13:13 PDT 2011
Author: becker
Date: 2011-05-15 15:13:13 -0700 (Sun, 15 May 2011)
New Revision: 18376
Modified:
mc/1D/hc/trunk/anisotropic_viscosity/Anisotropic_viscosity.c
mc/1D/hc/trunk/ggrd_grdtrack_util.c
Log:
Fixed typo in CitcomS version of anisotropic viscosity layered init.
Modified: mc/1D/hc/trunk/anisotropic_viscosity/Anisotropic_viscosity.c
===================================================================
--- mc/1D/hc/trunk/anisotropic_viscosity/Anisotropic_viscosity.c 2011-05-15 03:43:09 UTC (rev 18375)
+++ mc/1D/hc/trunk/anisotropic_viscosity/Anisotropic_viscosity.c 2011-05-15 22:13:13 UTC (rev 18376)
@@ -378,19 +378,25 @@
ggrd_read_anivisc_from_file(E);
#endif
break;
- case 6: /* tapered within layer */
+ case 6:
+ /*
+ tapered within layer
+ */
if(E->parallel.me == 0)
fprintf(stderr,"set_anisotropic_viscosity_at_element_level: setting orthotropic tapered, vis2 min %g\n",
E->viscosity.ani_vis2_factor);
- if(E->viscosity.anivisc_layer >= 0)myerror_s("set_anisotropic_viscosity_at_element_level: need to select layer",E);
+ if(E->viscosity.anivisc_layer >= 0)
+ myerror_s("set_anisotropic_viscosity_at_element_level: need to select layer",E);
ani_layer = -E->viscosity.anivisc_layer;
-#ifdef CitcomS_global_defs_h /* CitcomS */
- z_bottom = E->viscosity.zbase_layer[ani_layer-1];
+#ifdef CitcomS_global_defs_h
+ /* CitcomS */
+ z_bottom = E->sphere.ro-E->viscosity.zbase_layer[ani_layer-1];
if(ani_layer == 1)
z_top = E->sphere.ro;
else
- z_top = E->viscosity.zbase_layer[ani_layer-2];
-#else /* CU */
+ z_top = E->sphere.ro - E->viscosity.zbase_layer[ani_layer-2];
+#else
+ /* CU */
z_bottom = E->viscosity.zbase_layer[ani_layer-1];
if(ani_layer == 1)
z_top = E->segment.zzlayer[E->segment.zlayers-1];
@@ -405,9 +411,9 @@
elxlz = elx * elz;
for (j=1;j <= elz;j++){
#ifdef CitcomS_global_defs_h /* CitcomS */
- if(E->mat[m][j] == -E->viscosity.anivisc_layer){
+ if(E->mat[m][j] == ani_layer){
#else
- if(E->mat[j] == -E->viscosity.anivisc_layer){
+ if(E->mat[j] == ani_layer){
#endif
for(u=0.,inode=1;inode <= ends;inode++){ /* mean vertical coordinate */
#ifdef CitcomS_global_defs_h /* CitcomS */
@@ -422,16 +428,21 @@
#endif
}
u /= ends;
- /* do a log scale decrease of vis2 to ani_vis2_factor from bottom to top of layer */
+ /*
+ do a log scale decrease of vis2 to ani_vis2_factor from bottom to top of layer
+ */
vis2 = exp(log(E->viscosity.ani_vis2_factor) * (u-z_bottom)/(z_top-z_bottom));
//fprintf(stderr,"z %g (%g/%g) vis2 %g vis2_o %g frac %g\n",u,z_top,z_bottom,vis2, E->viscosity.ani_vis2_factor,(u-z_bottom)/(z_top-z_bottom));
- /* 1-eta_s/eta */
+ /*
+ 1-eta_s/eta
+ */
vis2 = 1 - vis2;
for (k=1;k <= ely;k++){
for (l=1;l <= elx;l++) {
/* eq.(1) */
el = j + (l-1) * elz + (k-1)*elxlz;
-#ifdef CitcomS_global_defs_h /* CitcomS */
+#ifdef CitcomS_global_defs_h
+ /* CitcomS */
xloc[0] = xloc[1] = xloc[2] = 0.0;
for(inode=1;inode <= ends;inode++){
off = E->ien[m][el].node[inode];
@@ -439,15 +450,21 @@
xloc[0] += rout[0];xloc[1] += rout[1];xloc[2] += rout[2];
}
xloc[0]/=ends;xloc[1]/=ends;xloc[2]/=ends;xyz2rtp(xloc[0],xloc[1],xloc[2],rout);
- calc_cbase_at_tp(rout[1],rout[2],base);convert_pvec_to_cvec(1.,0.,0.,base,rout);
+ /* r,t,p=(1,0,0) convert to Caretsian */
+ calc_cbase_at_tp(rout[1],rout[2],base);
+ convert_pvec_to_cvec(1.,0.,0.,base,rout);
n[0]=rout[0];n[1]=rout[1];n[2]=rout[2];
for(p=1;p <= vpts;p++){ /* assign to all integration points */
off = (el-1)*vpts + p;
E->EVI2[i][m][off] = vis2;
E->EVIn1[i][m][off] = n[0]; E->EVIn2[i][m][off] = n[1];E->EVIn3[i][m][off] = n[2];
+ //fprintf(stderr,"%g %g %g %g\n",n[0],n[1],n[2],vis2);
E->avmode[i][m][off] = CITCOM_ANIVISC_ORTHO_MODE;
}
-#else /* CU */
+#else
+ /*
+ CU
+ */
if(E->control.Rsphere){ /* director in r direction */
xloc[0] = xloc[1] = xloc[2] = 0.0;
for(inode=1;inode <= ends;inode++){
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2011-05-15 03:43:09 UTC (rev 18375)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2011-05-15 22:13:13 UTC (rev 18376)
@@ -104,7 +104,8 @@
g->east = g->west = g->south = g->north = 0.0;
g->is_three = is_three;
g->init = FALSE;
- if(ggrd_grdtrack_init(&g->west,&g->east,&g->south,&g->north,&g->f,&g->mm,
+ if(ggrd_grdtrack_init(&g->west,&g->east,
+ &g->south,&g->north,&g->f,&g->mm,
grdfile,&g->grd,&g->edgeinfo,
gmt_edgeinfo_string,&g->geographic_in,
pad,is_three,depth_file,&g->z,&g->nz,
@@ -120,7 +121,7 @@
for(i=0;i < g->nz;i++){ /* loop through layers */
//g->fmaxlim[i] = g->grd[i].z_min;
g->fmaxlim[i] = 0.0;
- for(j=0;j<g->mm;j++){ /* loop trough entries */
+ for(j=0;j < g->mm;j++){ /* loop trough entries */
tmp = fabs(g->f[i*g->mm+j]);
//if((g->f[i*g->mm+j] < g->bandlim) &&(g->f[i*g->mm+j] > g->fmaxlim[i]))
if((tmp < g->bandlim) && tmp > g->fmaxlim[i])
@@ -443,12 +444,17 @@
levels to go from depth (>0) to z (<0) */
struct GMT_BCR *loc_bcr)
#else
-int ggrd_grdtrack_init(double *west, double *east,double *south, double *north,
- float **f,int *mm,char *grdfile,struct GRD_HEADER **grd,
- struct GMT_EDGEINFO **edgeinfo,char *edgeinfo_string,
- ggrd_boolean *geographic_in,int *pad,ggrd_boolean three_d,
+int ggrd_grdtrack_init(double *west, double *east,
+ double *south, double *north,
+ float **f,int *mm,char *grdfile,
+ struct GRD_HEADER **grd,
+ struct GMT_EDGEINFO **edgeinfo,
+ char *edgeinfo_string,
+ ggrd_boolean *geographic_in,
+ int *pad,ggrd_boolean three_d,
char *dfile, float **z,int *nz,
- ggrd_boolean interpolant,ggrd_boolean verbose,
+ ggrd_boolean interpolant,
+ ggrd_boolean verbose,
ggrd_boolean change_depth_sign,
struct BCR *loc_bcr)
#endif
@@ -750,7 +756,8 @@
#endif
}
/* Set boundary conditions */
- GMT_boundcond_set ((*grd+i), (*edgeinfo+i), pad, (*f+i*(*mm)));
+ GMT_boundcond_set ((*grd+i), (*edgeinfo+i), pad,
+ (*f+i*(*mm)));
} /* end layer loop */
if(verbose){
ggrd_print_layer_avg(*f,*z,mx,my,*nz,stderr,pad);
@@ -758,7 +765,9 @@
return 0;
}
-void ggrd_print_layer_avg(float *x,float *z,int nx, int ny, int m,FILE *out,
+
+void ggrd_print_layer_avg(float *x,float *z,int nx, int ny,
+ int m,FILE *out,
GMT_LONG *pad) /* >= 4.5.1 */
//int *pad)
{
@@ -768,7 +777,8 @@
if(pad[0]+pad[1]+pad[2]+pad[3] == 0){
for(i=0;i < m;i++){
fprintf(stderr,"ggrd_grdtrack_init: layer %3i at depth %11g, mean: %11g rms: %11g\n",
- i+1,z[i],ggrd_gt_mean((x+i*nxny),nxny),ggrd_gt_rms((x+i*nxny),nxny));
+ i+1,z[i],ggrd_gt_mean((x+i*nxny),nxny),
+ ggrd_gt_rms((x+i*nxny),nxny));
}
}else{
@@ -782,7 +792,8 @@
for(k=pad[0];k < xl;k++,l++)
tmp[l] = x[i*nxny + j * nx + k];
fprintf(stderr,"ggrd_grdtrack_init: layer %3i at depth %11g, mean: %11g rms: %11g\n",
- i+1,z[i],ggrd_gt_mean(tmp,nxnyr),ggrd_gt_rms(tmp,nxnyr));
+ i+1,z[i],ggrd_gt_mean(tmp,nxnyr),
+ ggrd_gt_rms(tmp,nxnyr));
}
free(tmp);
}
More information about the CIG-COMMITS
mailing list