[cig-commits] r17292 - mc/3D/CitcomCU/trunk/src
becker at geodynamics.org
becker at geodynamics.org
Sun Oct 17 19:07:28 PDT 2010
Author: becker
Date: 2010-10-17 19:07:28 -0700 (Sun, 17 Oct 2010)
New Revision: 17292
Modified:
mc/3D/CitcomCU/trunk/src/Anisotropic_viscosity.c
mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
Log:
beta version of anisotropic viscosity working.
Modified: mc/3D/CitcomCU/trunk/src/Anisotropic_viscosity.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Anisotropic_viscosity.c 2010-10-18 00:58:09 UTC (rev 17291)
+++ mc/3D/CitcomCU/trunk/src/Anisotropic_viscosity.c 2010-10-18 02:07:28 UTC (rev 17292)
@@ -229,10 +229,11 @@
}
void set_anisotropic_viscosity_at_element_level(struct All_variables *E, int init_stage)
{
- int i,k,l,off,nel;
- double vis2,n[3],u,v,s,r;
+ int i,j,k,l,m,off,nel,elx,ely,elz,inode,elxlz,el,ani_layer;
+ double vis2,n[3],u,v,s,r,xloc[3],z_top,z_bottom;
+ float base[9],rout[3];
const int vpts = vpoints[E->mesh.nsd];
-
+ const int ends = enodes[E->mesh.nsd];
if(init_stage){
if(E->parallel.me == 0)
fprintf(stderr,"set_anisotropic_viscosity: allowing for %s viscosity\n",
@@ -300,6 +301,76 @@
#endif
ggrd_read_anivisc_from_file(E);
break;
+ 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("set_anisotropic_viscosity_at_element_level: need to select layer",E);
+ ani_layer = -E->viscosity.anivisc_layer;
+ z_bottom = E->viscosity.zbase_layer[ani_layer-1];
+ if(ani_layer == 1)
+ z_top = E->segment.zzlayer[E->segment.zlayers-1];
+ else
+ z_top = E->viscosity.zbase_layer[ani_layer-2];
+
+ for(i=E->mesh.levmin;i <= E->mesh.levmax;i++){
+ elx = E->lmesh.ELX[i];elz = E->lmesh.ELZ[i];ely = E->lmesh.ELY[i];
+ elxlz = elx * elz;
+ for (j=1;j <= elz;j++)
+ if(E->mat[j] == -E->viscosity.anivisc_layer){
+
+ for(u=0.,inode=1;inode <= ends;inode++){ /* mean vertical coordinate */
+ off = E->ien[j].node[inode];
+ if(E->control.Rsphere)
+ u += E->SX[3][off];
+ else
+ u += E->X[3][off];
+ }
+ u /= ends;
+ /* 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 */
+ 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;
+ if(E->control.Rsphere){ /* director in r direction */
+ xloc[0] = xloc[1] = xloc[2] = 0.0;
+ for(inode=1;inode <= ends;inode++){
+ off = E->ien[el].node[inode];
+ rtp2xyz((float)E->SX[3][off],(float)E->SX[1][off],(float)E->SX[2][off],rout);
+ 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);
+ n[0]=rout[0];n[1]=rout[1];n[2]=rout[2];
+ }else{ /* director in z direction */
+ n[0] = 0.;n[1] = 0.;n[2] = 1.;
+ }
+ for(m=1;m <= vpts;m++){ /* assign to all integration points */
+ off = (el-1)*vpts + m;
+ E->EVI2[i][off] = vis2;
+ E->EVIn1[i][off] = n[0]; E->EVIn2[i][off] = n[1];E->EVIn3[i][off] = n[2];
+ E->avmode[i][off] = CITCOM_ANIVISC_ORTHO_MODE;
+ }
+ }
+ }
+ }else{ /* outside layer = isotropic */
+ for (k=1;k <= ely;k++)
+ for (l=1;l <= elx;l++){
+ el = j + (l-1) * elz + (k-1)*elxlz;
+ for(m=1;m <= vpts;m++){ /* assign to all integration points */
+ off = (el-1)*vpts + m;
+ E->EVI2[i][off] = 0;
+ E->EVIn1[i][off] = 1; E->EVIn2[i][off] = 0;E->EVIn3[i][off] = 0;
+ E->avmode[i][off] = CITCOM_ANIVISC_ORTHO_MODE;
+ }
+ }
+ }
+ } /* end multigrid level loop */
+ break;
default:
fprintf(stderr,"set_anisotropic_viscosity_at_element_level: anisotropic_init %i undefined\n",
E->viscosity.anisotropic_init);
@@ -327,7 +398,8 @@
if((E->monitor.solution_cycles > 0) || (E->monitor.visc_iter_count > 0))
align_director_with_ISA_for_element(E,CITCOM_ANIVISC_MIXED_ALIGN);
break;
- default: /* default, no action */
+ default: /* default, no further modification of
+ anisotropy */
break;
}
}
Modified: mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2010-10-18 00:58:09 UTC (rev 17291)
+++ mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2010-10-18 02:07:28 UTC (rev 17292)
@@ -578,7 +578,7 @@
MPI_Status mpi_stat;
int mpi_rc;
int mpi_inmsg, mpi_success_message = 1;
- int el,i,j,k,l,inode,i1,i2,elxlz,elxlylz,ind,nel;
+ int el,i,j,k,l,inode,i1,i2,elxlz,ind,nel;
int llayer,nox,noy,noz,level,lselect,idim,elx,ely,elz;
char gmt_string[10],char_dummy;
double vis2,log_vis,ntheta,nphi,nr;
@@ -591,15 +591,13 @@
const int ends = enodes[dims];
FILE *in;
struct ggrd_gt *vis2_grd,*ntheta_grd,*nphi_grd,*nr_grd;
- const int vpts = vpoints[E->mesh.nsd];
+ const int vpts = vpoints[dims];
static int init = FALSE;
-
+
nox=E->mesh.nox;noy=E->mesh.noy;noz=E->mesh.noz;
elx=E->lmesh.elx;elz=E->lmesh.elz;ely=E->lmesh.ely;
elxlz = elx * elz;
- elxlylz = elxlz * ely;
-
if(E->viscosity.allow_anisotropic_viscosity == 0)
myerror("ggrd_read_anivisc_from_file: called, but allow_anisotropic_viscosity is FALSE?!",E);
if(init)
@@ -712,13 +710,13 @@
*/
if(E->control.CART3D){
rout[0]=rout[1]=rout[2]=0.0;
- for(inode=1;inode <= 4;inode++){
+ for(inode=1;inode <= ends;inode++){
ind = E->ien[el].node[inode];
rout[0] += E->X[1][ind];
rout[1] += E->X[2][ind];
rout[2] += E->X[3][ind];
}
- rout[0]/=4.;rout[1]/=4.;rout[2]/=4.;
+ rout[0]/=ends;rout[1]/=ends;rout[2]/=ends;
if(!ggrd_grdtrack_interpolate_xy((double)rout[0],(double)rout[1],vis2_grd,&log_vis,FALSE)){
fprintf(stderr,"ggrd_read_anivisc_from_file: interpolation error at x: %g y: %g\n",
rout[0],rout[1]);
@@ -744,7 +742,7 @@
}else{
/* spherical */
xloc[0] = xloc[1] = xloc[2] = 0.0;
- for(inode=1;inode <= 4;inode++){
+ for(inode=1;inode <= ends;inode++){
ind = E->ien[el].node[inode];
rtp2xyz((float)E->SX[3][ind],(float)E->SX[1][ind],
(float)E->SX[2][ind],rout);
@@ -752,7 +750,7 @@
xloc[1] += rout[1];
xloc[2] += rout[2];
}
- xloc[0]/=4.;xloc[1]/=4.;xloc[2]/=4.;
+ xloc[0]/=ends;xloc[1]/=ends;xloc[2]/=ends;
xyz2rtp(xloc[0],xloc[1],xloc[2],rout); /* convert to spherical */
/* vis2 */
Modified: mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2010-10-18 00:58:09 UTC (rev 17291)
+++ mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2010-10-18 02:07:28 UTC (rev 17292)
@@ -190,9 +190,10 @@
1: random
2: read in director orientation
and log10(eta_s/eta)
- 3: align with velocity
- 4: align with ISA
- 5: align mixed depending on deformation state
+ 3: align with velocity, use ani_vis2_factor for eta_s/eta
+ 4: align with ISA, use ani_vis2_factor for eta_s/eta
+ 5: align mixed depending on deformation state, use ani_vis2_factor for eta_s/eta
+ 6: use radially aligned director and taper eta_s/eta from base (1) to top of layer (ani_vis2_factor)
*/
input_string("anisotropic_init_dir",(E->viscosity.anisotropic_init_dir),"",m); /* directory
@@ -203,6 +204,8 @@
input_int("anivisc_layer",&(E->viscosity.anivisc_layer),"1",m); /* >0: assign to layers on top of anivisc_layer
<0: assign to layer = anivisc_layer
*/
+ if((E->viscosity.anisotropic_init == 6) && (E->viscosity.anivisc_layer >= 0))
+ myerror("anisotropic init mode 6 requires selection of layer where anisotropy applies",E);
input_boolean("anivisc_start_from_iso",
&(E->viscosity.anivisc_start_from_iso),"on",m); /* start
@@ -215,7 +218,7 @@
E->viscosity.anivisc_start_from_iso = TRUE;
}
/* ratio between weak and strong direction */
- input_double("ani_vis2_factor",&(E->viscosity.ani_vis2_factor),"0.1",m);
+ input_double("ani_vis2_factor",&(E->viscosity.ani_vis2_factor),"1.0",m);
}
More information about the CIG-COMMITS
mailing list