[cig-commits] r17174 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Mon Sep 6 18:40:03 PDT 2010
Author: becker
Date: 2010-09-06 18:40:02 -0700 (Mon, 06 Sep 2010)
New Revision: 17174
Modified:
mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
mc/3D/CitcomS/trunk/lib/Nodal_mesh.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
Log:
Checking in current version because of need for debugging.
Revision 17172, compared to (presumably) 17153 give strange multigrid
convergence behavior. Checking why that is now.
Modified: mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c 2010-09-06 16:21:27 UTC (rev 17173)
+++ mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c 2010-09-07 01:40:02 UTC (rev 17174)
@@ -55,7 +55,8 @@
output: D[0,...,5][0,...,5] constitutive matrix
- input: delta_vis difference in viscosity from isotropic viscosity, which is set to unity
+
+ input: delta_vis difference in viscosity from isotropic viscosity (set to unity here)
n[0,..,2]: director orientation, in cartesian
@@ -343,7 +344,8 @@
for(j2=0;j2<3;j2++)
for(j3=0;j3<3;j3++)
for(j4=0;j4<3;j4++)
- c4c[i1][i2][i3][i4] += r[i1][j1] * r[i2][j2]* r[i3][j3]* r[i4][j4] * c4[j1][j2][j3][j4];
+ c4c[i1][i2][i3][i4] += r[i1][j1] * r[i2][j2] *
+ r[i3][j3]* r[i4][j4] * c4[j1][j2][j3][j4];
}
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-09-06 16:21:27 UTC (rev 17173)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-09-07 01:40:02 UTC (rev 17174)
@@ -1318,7 +1318,7 @@
read in anisotropic viscosity from a directory which holds
-vis2.grd for the viscosity factors
+vis2.grd for the viscosity factors (1 - eta_S/eta)
nr.grd, nt.grd, np.grd for the directors
@@ -1339,7 +1339,6 @@
const int ends = enodes[dims];
FILE *in;
struct ggrd_gt *vis2_grd,*ntheta_grd,*nphi_grd,*nr_grd;
- const int init_layer = 1; /* initialize all elements in top layer */
const int vpts = vpoints[E->mesh.nsd];
@@ -1395,9 +1394,16 @@
if(E->parallel.me==0)
- fprintf(stderr,"ggrd_read_anivisc_from_file: initializing, assigning to all above %g km, input is %s\n",
- E->data.radius_km*E->viscosity.zbase_layer[E->control.ggrd.mat_control-1],
- (is_global)?("global"):("regional"));
+ if(E->viscosity.anivisc_layer > 0)
+ fprintf(stderr,"ggrd_read_anivisc_from_file: initializing, assigning to all elements above %g km, input is %s\n",
+ E->data.radius_km*E->viscosity.zbase_layer[E->viscosity.anivisc_layer - 1],
+ (is_global)?("global"):("regional"));
+ else
+ fprintf(stderr,"ggrd_read_anivisc_from_file: initializing, assigning to all elements between %g and %g km, input is %s\n",
+ E->data.radius_km*((E->viscosity.anivisc_layer<-1)?(E->viscosity.zbase_layer[-E->viscosity.anivisc_layer - 2]):(0)),
+ E->data.radius_km*E->viscosity.zbase_layer[-E->viscosity.anivisc_layer - 1],
+ (is_global)?("global"):("regional"));
+
/*
read viscosity ratio, and east/north direction of normal azimuth
*/
@@ -1437,7 +1443,8 @@
*/
for (m=1;m <= E->sphere.caps_per_proc;m++) {
for (j=1;j <= elz;j++) { /* this assumes a regular grid sorted as in (1)!!! */
- if(E->mat[m][j] <= init_layer){
+ if(((E->viscosity.anivisc_layer > 0)&&(E->mat[m][j] <= E->viscosity.anivisc_layer))||
+ ((E->viscosity.anivisc_layer < 0)&&(E->mat[m][j] == -E->viscosity.anivisc_layer))){
/* within top layers */
for (k=1;k <= ely;k++){
for (i=1;i <= elx;i++) {
@@ -1449,7 +1456,9 @@
xloc[1] = xloc[2] = xloc[3] = 0.0;
for(inode=1;inode <= 4;inode++){
ind = E->ien[m][el].node[inode];
- xloc[1] += E->x[m][1][ind];xloc[2] += E->x[m][2][ind];xloc[3] += E->x[m][3][ind];
+ xloc[1] += E->x[m][1][ind];
+ xloc[2] += E->x[m][2][ind];
+ xloc[3] += E->x[m][3][ind];
}
xloc[1]/=4.;xloc[2]/=4.;xloc[3]/=4.;
xyz2rtpd(xloc[1],xloc[2],xloc[3],rout); /* convert to spherical */
@@ -1482,16 +1491,20 @@
rout[2]*180/M_PI,90-rout[1]*180/M_PI);
parallel_process_termination();
}
- /*
- convert from n_east,n_north to Cartesian vector,
- assuming the director is in the horizontal,
- i.e. n_r = 0
- */
- nlen = sqrt(nphi*nphi+ntheta*ntheta+nr*nr); /* correct, because
- interpolation might have
- screwed up initialization */
+ nlen = sqrt(nphi*nphi + ntheta*ntheta + nr*nr); /* correct,
+ because
+ interpolation
+ might
+ have
+ screwed
+ up
+ initialization */
nphi /= nlen; ntheta /= nlen;nr /= nlen;
- calc_cbase_at_tp_d(rout[1],rout[2],base);
+ calc_cbase_at_tp_d(rout[1],rout[2],base); /* convert from
+ spherical
+ coordinates
+ to
+ Cartesian */
convert_pvec_to_cvec_d(nr,ntheta,nphi,base,cvec);
for(l=1;l <= vpts;l++){ /* assign to all integration points */
ind = (el-1)*vpts + l;
Modified: mc/3D/CitcomS/trunk/lib/Nodal_mesh.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Nodal_mesh.c 2010-09-06 16:21:27 UTC (rev 17173)
+++ mc/3D/CitcomS/trunk/lib/Nodal_mesh.c 2010-09-07 01:40:02 UTC (rev 17174)
@@ -245,7 +245,7 @@
const int nsd=E->mesh.nsd;
const int vpts=vpoints[nsd];
const int ends=enodes[nsd];
- double temp_visc, weight;
+ double temp_visc;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.NNO[lev];i++)
@@ -352,7 +352,6 @@
const int ends=enodes[nsd];
double temp_visc;
-
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(e=1;e<=E->lmesh.NEL[lev];e++)
for(i=1;i<=vpts;i++) {
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2010-09-06 16:21:27 UTC (rev 17173)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2010-09-07 01:40:02 UTC (rev 17174)
@@ -82,6 +82,9 @@
ggrd
type
init */
+ 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
+ */
}
Modified: mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2010-09-06 16:21:27 UTC (rev 17173)
+++ mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2010-09-07 01:40:02 UTC (rev 17174)
@@ -39,6 +39,7 @@
#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
int anisotropic_init; /* 0: isotropic, 1: random, 2: from file */
char anisotropic_init_dir[1000];
+ int anivisc_layer; /* layer to assign anisotropic viscosity to for mode 2 */
#endif
char STRUCTURE[20]; /* which option to determine viscosity field, one of .... */
More information about the CIG-COMMITS
mailing list