[cig-commits] r17185 - in mc/3D/CitcomS/trunk: . lib
becker at geodynamics.org
becker at geodynamics.org
Mon Sep 13 00:13:07 PDT 2010
Author: becker
Date: 2010-09-13 00:13:07 -0700 (Mon, 13 Sep 2010)
New Revision: 17185
Modified:
mc/3D/CitcomS/trunk/configure.ac
mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
mc/3D/CitcomS/trunk/lib/Element_calculations.c
mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Output_gzdir.c
mc/3D/CitcomS/trunk/lib/Solver_multigrid.c
mc/3D/CitcomS/trunk/lib/Topo_gravity.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
Log:
Made sure current version compiled nicely without additional flags.
Anisotropic viscosity is set up more flexibly.
Modified: mc/3D/CitcomS/trunk/configure.ac
===================================================================
--- mc/3D/CitcomS/trunk/configure.ac 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/configure.ac 2010-09-13 07:13:07 UTC (rev 17185)
@@ -200,8 +200,9 @@
if test -n "$GMTHOME"; then
LDFLAGS="$LDFLAGS -L$GMTHOME/lib"
fi
- AC_SEARCH_LIBS([GMT_grd_init], [gmt], [], [AC_MSG_ERROR([GMT library not found, required by ggrd. Try setting environment variable GMTHOME.])])
+ AC_SEARCH_LIBS([GMT_grd_init], [gmt], [], [AC_MSG_ERROR([GMT library not found, required by ggrd. Try setting environment variable GMTHOME.])], [-lm])
+
# Checking HC ggrd library
if test -n "$HC_HOME"; then
if test -n "$ARCH"; then
@@ -210,7 +211,7 @@
LDFLAGS="$LDFLAGS -L$HC_HOME/objects"
fi
fi
- AC_SEARCH_LIBS([ggrd_init_master], [ggrd], [], [AC_MSG_ERROR([HC ggrd library not found. Try setting environment variable HC_HOME.])])
+ AC_SEARCH_LIBS([ggrd_init_master], [ggrd], [], [AC_MSG_ERROR([HC ggrd library not found. Try setting environment variable HC_HOME.])],[-lm])
fi
AC_SEARCH_LIBS([sqrt], [m])
Modified: mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c 2010-09-13 07:13:07 UTC (rev 17185)
@@ -33,11 +33,12 @@
anisotropic viscosity following Muehlhaus, Moresi, Hobbs and Dufour
(PAGEOPH, 159, 2311, 2002)
+ tranverse isotropy following Han and Wahr (PEPI, 102, 33, 1997)
*/
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
#include <math.h>
#include "element_definitions.h"
@@ -48,17 +49,91 @@
void calc_cbase_at_tp_d(double , double , double *);
#define CITCOM_DELTA(i,j) ((i==j)?(1.0):(0.0))
+
+
/*
+transversely isotropic viscosity following Han and Wahr
-compute a cartesian anisotropic viscosity matrix
+\nu_1 = isotropic viscosity, applies for e_31, e_23
+\nu_2 = weak anisotropy, applies for e_31, e_32
+\eta_1 = normal viscosity, (\eta_1+2\nu_1) control e_11, e_22
+\eta_2 = normal viscosity, (\eta_2+2\nu_2) = 2\eta_1 + 2\nu_1, controls e_33
+we use (for consistency with anisotropic viscosity)
+
+Delta = 1-\nu_2/\nu_1
+
+and
+
+\Gamma, such that \eta_1 = \Gamma \nu_1
+
+\nu_1 is the reference, isotropic viscosity, set to unity here, i.e.
+
+\nu_2 = 1 - \Delta ; \eta_1 = \Gamma ; (\eta_2 = 2 (\Gamma-\Delta)); for isotropy \Delta = 0, \Gamma = 0
+
+n[3] is the cartesian direction into which the weak shear points
+(ie. routine will rotate the 3 axis into the n diretion)
+
+*/
+void get_constitutive_ti_viscosity(double D[6][6], double delta_vis, double gamma_vis,
+ double n[3], int convert_to_spherical,
+ double theta, double phi)
+{
+ double nlen,delta_vis2;
+ int ani;
+ /*
+ make sure things are normalized (n[3] might come from interpolation)
+ */
+ nlen = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
+ if((nlen < 0.95)||(nlen > 1.05)){
+ fprintf(stderr,"get_constitutive_ti_viscosity: error: nlen: %g\n",nlen);
+ parallel_process_termination();
+ }
+
+ /* isotropic part, in units of iso_visc */
+ zero_6x6(D);
+ D[0][0] = 2.0; /* xx = tt*/
+ D[1][1] = 2.0; /* yy = pp */
+ D[2][2] = 2.0; /* zz = rr */
+ D[3][3] = 1.0; /* xy = tp */
+ D[4][4] = 1.0; /* xz = rt */
+ D[5][5] = 1.0; /* yz = rp */
+
+ ani = FALSE;
+ if((fabs(delta_vis) > 3e-15) || (fabs(gamma_vis) > 3e-15)){
+ ani = TRUE;
+ /* get Cartesian anisotropy correction matrix */
+ D[0][0] += gamma_vis;
+ D[1][0] = D[0][1] = D[1][1] = D[0][0];
+ D[2][2] += 2.*gamma_vis;
+ D[4][4] -= delta_vis;
+ D[5][5] = D[4][4];
+ //print_6x6_mat(stderr,D);
+ rotate_ti6x6_to_director(D,n); /* rotate, can use same mat for 6x6 */
+ //print_6x6_mat(stderr,D);fprintf(stderr,"\n");
+ }
+ if(ani && convert_to_spherical){
+ conv_cart6x6_to_spherical(D,theta,phi,D); /* rotate, can use in place */
+ }
+}
+
+/*
+
+
+compute a cartesian orthotropic anisotropic viscosity matrix (and
+rotate it into CitcomS spherical, if requested)
+
+viscosity is characterized by a eta_S (weak) viscosity in a shear
+plane, to which the director is normal
+
+
output: D[0,...,5][0,...,5] constitutive matrix
input: delta_vis difference in viscosity from isotropic viscosity (set to unity here)
- n[0,..,2]: director orientation, in cartesian
+ n[0,..,2]: director orientation, specify in cartesian
where delta_vis = (1 - eta_S/eta)
@@ -70,8 +145,8 @@
double theta, double phi)
{
double nlen,delta_vis2;
- double delta[3][3][3][3],deltac[3][3][3][3];
-
+ double delta[3][3][3][3];
+ int ani;
/*
make sure things are normalized (n[3] might come from interpolation)
*/
@@ -80,16 +155,9 @@
fprintf(stderr,"get_constitutive_orthotropic_viscosity: error: nlen: %g\n",nlen);
parallel_process_termination();
}
-
- /* zero out most of D matrix */
- D[0][1] = D[0][2] = D[0][3] = D[0][4] = D[0][5] = 0.;
- D[1][0] = D[1][2] = D[1][3] = D[1][4] = D[1][5] = 0.;
- D[2][0] = D[2][1] = D[2][3] = D[2][4] = D[2][5] = 0.;
- D[3][0] = D[3][1] = D[3][2] = D[3][4] = D[3][5] = 0.;
- D[4][0] = D[4][1] = D[4][2] = D[4][3] = D[4][5] = 0.;
- D[5][0] = D[5][1] = D[5][2] = D[5][3] = D[5][4] = 0.;
-
+ ani=FALSE;
/* isotropic part, in units of iso_visc */
+ zero_6x6(D);
D[0][0] = 2.0; /* xx = tt*/
D[1][1] = 2.0; /* yy = pp */
D[2][2] = 2.0; /* zz = rr */
@@ -97,20 +165,11 @@
D[4][4] = 1.0; /* xz = rt */
D[5][5] = 1.0; /* yz = rp */
-
- if(fabs(delta_vis) > 5e-15){
- /* get Cartesian anisotropy matrix */
- if(convert_to_spherical){
- get_delta(deltac,n); /* get anisotropy tensor, \Delta of
+ /* get Cartesian anisotropy matrix */
+ if(fabs(delta_vis) > 3e-15){
+ ani = TRUE;
+ get_orth_delta(delta,n); /* get anisotropy tensor, \Delta of
Muehlhaus et al. (2002) */
- conv_cart4x4_to_spherical(deltac,theta,phi,delta); /* rotate
- into
- CitcomS
- spherical
- system */
- }else{
- get_delta(delta,n);
- }
delta_vis2 = 2.0*delta_vis;
/* s_xx = s_tt */
D[0][0] -= delta_vis2 * delta[0][0][0][0]; /* * e_xx */
@@ -155,8 +214,11 @@
D[5][4] -= delta_vis * (delta[1][2][0][2]+delta[1][2][2][0]);
D[5][5] -= delta_vis * (delta[1][2][1][2]+delta[1][2][2][1]);
}
-
-
+ if(ani && convert_to_spherical){
+ //print_6x6_mat(stderr,D);
+ conv_cart6x6_to_spherical(D,theta,phi,D); /* rotate, can use same mat for 6x6 */
+ //print_6x6_mat(stderr,D);fprintf(stderr,"\n");
+ }
}
void set_anisotropic_viscosity_at_element_level(struct All_variables *E, int init_stage)
@@ -165,9 +227,9 @@
double vis2,n[3],u,v,s,r;
const int vpts = vpoints[E->mesh.nsd];
- if(E->viscosity.allow_orthotropic_viscosity){
+ if(E->viscosity.allow_anisotropic_viscosity){
if(init_stage){
- if(E->viscosity.orthotropic_viscosity_init)
+ if(E->viscosity.anisotropic_viscosity_init)
myerror(E,"anisotropic viscosity should not be initialized twice?!");
/* first call */
/* initialize anisotropic viscosity at element level, nodes will
@@ -236,7 +298,7 @@
parallel_process_termination();
break;
}
- E->viscosity.orthotropic_viscosity_init = TRUE;
+ E->viscosity.anisotropic_viscosity_init = TRUE;
/* end initialization stage */
}else{
/* standard operation every time step */
@@ -245,8 +307,8 @@
}
} /* end anisotropic viscosity branch */
}
-#endif
+
void normalize_director_at_nodes(struct All_variables *E,float **n1,float **n2, float **n3, int lev)
{
int n,m;
@@ -276,8 +338,11 @@
}
}
/*
+
+convert cartesian fourth order tensor (input c) to spherical, CitcomS
+format (output p)
-convert cartesian voigt matrix to spherical, CitcomS format
+c and p cannot be the same matrix
1: t 2: p 3: r
@@ -287,25 +352,61 @@
void conv_cart4x4_to_spherical(double c[3][3][3][3], double theta, double phi, double p[3][3][3][3])
{
- double rot[3][3],base[9];
+ double rot[3][3];
+ get_citcom_spherical_rot(theta,phi,rot);
+ rot_4x4(c,rot,p);
+}
+
+/* convert [6][6] (input c) in cartesian to citcom spherical (output
+ p)
+
+ c and p can be the same amtrix
+
+*/
+void conv_cart6x6_to_spherical(double c[6][6], double theta, double phi, double p[6][6])
+{
+ double c4[3][3][3][3],p4[3][3][3][3],rot[3][3];
+ get_citcom_spherical_rot(theta,phi,rot);
+ c4fromc6(c4,c);
+ rot_4x4(c4,rot,p4);
+ c6fromc4(p,p4);
+}
+/*
+
+rotate 6x6 D matrix with preferred axis aligned with z to the
+Cartesian director orientation, in place
+
+ */
+void rotate_ti6x6_to_director(double D[6][6],double n[3])
+{
+ double a[3][3][3][3],b[3][3][3][3],rot[3][3];
+ rot[0][0] = rot[0][1] =
+ rot[1][0] = rot[1][1] =
+ rot[2][0] = rot[2][1] = 0.0;
+ rot[0][0] = n[0]; rot[1][0] = n[1]; rot[2][0] = n[2];
+ c4fromc6(a,D);
+ rot_4x4(a,rot,b);
+ c6fromc4(D,b);
+}
+
+void get_citcom_spherical_rot(double theta, double phi, double rot[3][3]){
+ double base[9];
calc_cbase_at_tp_d(theta,phi, base); /* compute cartesian basis at
theta, phi location */
rot[0][0] = base[3];rot[0][1] = base[4];rot[0][2] = base[5]; /* theta */
rot[1][0] = base[6];rot[1][1] = base[7];rot[1][2] = base[8]; /* phi */
rot[2][0] = base[0];rot[2][1] = base[1];rot[2][2] = base[2]; /* r */
- //fprintf(stderr,"%g %g ; %g %g %g ; %g %g %g ; %g %g %g\n\n",theta,phi,rot[0][0],rot[0][1],rot[0][2],rot[1][0],rot[1][1],rot[1][2],rot[2][0],rot[2][1],rot[2][2]);
- rot_4x4(c,rot,p);
- //if(E->parallel.me==0)print_6x6_mat(stderr,p);
- //if(E->parallel.me==0)fprintf(stderr,"\n\n");
+ //fprintf(stderr,"%g %g ; %g %g %g ; %g %g %g ; %g %g %g\n\n",
+ //theta,phi,rot[0][0],rot[0][1],rot[0][2],rot[1][0],rot[1][1],rot[1][2],rot[2][0],rot[2][1],rot[2][2]);
}
-
/*
-get anisotropy matrix following Muehlhaus
+get fourth order anisotropy tensor for orthotropic viscosity from
+Muehlhaus et al. (2002)
- */
-void get_delta(double d[3][3][3][3],double n[3])
+*/
+void get_orth_delta(double d[3][3][3][3],double n[3])
{
int i,j,k,l;
double tmp;
@@ -325,17 +426,18 @@
}
-/* rotate fourth order tensor */
+/*
+ rotate fourth order tensor
+ c4 and c4c cannot be the same matrix
+
+*/
void rot_4x4(double c4[3][3][3][3], double r[3][3], double c4c[3][3][3][3])
{
int i1,i2,i3,i4,j1,j2,j3,j4;
- for(i1=0;i1<3;i1++)
- for(i2=0;i2<3;i2++)
- for(i3=0;i3<3;i3++)
- for(i4=0;i4<3;i4++)
- c4c[i1][i2][i3][i4] = 0.0;
-
+
+ zero_4x4(c4c);
+
for(i1=0;i1<3;i1++)
for(i2=0;i2<3;i2++)
for(i3=0;i3<3;i3++)
@@ -348,7 +450,43 @@
r[i3][j3]* r[i4][j4] * c4[j1][j2][j3][j4];
}
+void zero_6x6(double a[6][6])
+{
+ int i,j;
+ for(i=0;i<6;i++)
+ for(j=0;j<6;j++)
+ a[i][j] = 0.;
+
+}
+void zero_4x4(double a[3][3][3][3])
+{
+ int i1,i2,i3,i4;
+ for(i1=0;i1<3;i1++)
+ for(i2=0;i2<3;i2++)
+ for(i3=0;i3<3;i3++)
+ for(i4=0;i4<3;i4++)
+ a[i1][i2][i3][i4] = 0.0;
+
+}
+void copy_4x4(double a[3][3][3][3], double b[3][3][3][3])
+{
+ int i1,i2,i3,i4;
+ for(i1=0;i1<3;i1++)
+ for(i2=0;i2<3;i2++)
+ for(i3=0;i3<3;i3++)
+ for(i4=0;i4<3;i4++)
+ b[i1][i2][i3][i4] = a[i1][i2][i3][i4];
+}
+void copy_6x6(double a[6][6], double b[6][6])
+{
+
+ int i1,i2;
+ for(i1=0;i1<6;i1++)
+ for(i2=0;i2<6;i2++)
+ b[i1][i2] = a[i1][i2];
+}
+
void print_6x6_mat(FILE *out, double c[6][6])
{
int i,j;
@@ -358,4 +496,104 @@
fprintf(out,"\n");
}
}
+/*
+ create a fourth order tensor representation from the voigt
+ notation, assuming only upper half is filled in
+ */
+void c4fromc6(double c4[3][3][3][3],double c[6][6])
+{
+ int i,j;
+
+ c4[0][0][0][0] = c[0][0];
+ c4[0][0][1][1] = c[0][1];
+ c4[0][0][2][2] = c[0][2];
+ c4[0][0][0][1] = c4[0][0][1][0] = c[0][3];
+ c4[0][0][0][2] = c4[0][0][2][0] = c[0][4];
+ c4[0][0][1][2] = c4[0][0][2][1] = c[0][5];
+
+ c4[1][1][0][0] = c[0][1];
+ c4[1][1][1][1] = c[1][1];
+ c4[1][1][2][2] = c[1][2];
+ c4[1][1][0][1] = c4[1][1][1][0] = c[1][3];
+ c4[1][1][0][2] = c4[1][1][2][0] = c[1][4];
+ c4[1][1][1][2] = c4[1][1][2][1] = c[1][5];
+
+ c4[2][2][0][0] = c[0][2];
+ c4[2][2][1][1] = c[1][2];
+ c4[2][2][2][2] = c[2][2];
+ c4[2][2][0][1] = c4[2][2][1][0] = c[2][3];
+ c4[2][2][0][2] = c4[2][2][2][0] = c[2][4];
+ c4[2][2][1][2] = c4[2][2][2][1] = c[2][5];
+
+ c4[0][1][0][0] = c[0][3];
+ c4[0][1][1][1] = c[1][3];
+ c4[0][1][2][2] = c[2][3];
+ c4[0][1][0][1] = c4[0][1][1][0] = c[3][3];
+ c4[0][1][0][2] = c4[0][1][2][0] = c[3][4];
+ c4[0][1][1][2] = c4[0][1][2][1] = c[3][5];
+
+ c4[0][2][0][0] = c[0][4];
+ c4[0][2][1][1] = c[1][4];
+ c4[0][2][2][2] = c[2][4];
+ c4[0][2][0][1] = c4[0][2][1][0] = c[3][4];
+ c4[0][2][0][2] = c4[0][2][2][0] = c[4][4];
+ c4[0][2][1][2] = c4[0][2][2][1] = c[4][5];
+
+ c4[1][2][0][0] = c[0][5];
+ c4[1][2][1][1] = c[1][5];
+ c4[1][2][2][2] = c[2][5];
+ c4[1][2][0][1] = c4[1][2][1][0] = c[3][5];
+ c4[1][2][0][2] = c4[1][2][2][0] = c[4][5];
+ c4[1][2][1][2] = c4[1][2][2][1] = c[5][5];
+
+ /* assign the symmetric diagonal terms */
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++){
+ c4[1][0][i][j] = c4[0][1][i][j];
+ c4[2][0][i][j] = c4[0][2][i][j];
+ c4[2][1][i][j] = c4[1][2][i][j];
+ }
+
+}
+void c6fromc4(double c[6][6],double c4[3][3][3][3])
+{
+ int i,j;
+
+ c[0][0] = c4[0][0][0][0];
+ c[0][1] = c4[0][0][1][1];
+ c[0][2] = c4[0][0][2][2];
+ c[0][3] = c4[0][0][0][1];
+ c[0][4] = c4[0][0][0][2];
+ c[0][5] = c4[0][0][1][2];
+
+ c[1][1] = c4[1][1][1][1];
+ c[1][2] = c4[1][1][2][2];
+ c[1][3] = c4[1][1][0][1];
+ c[1][4] = c4[1][1][0][2];
+ c[1][5] = c4[1][1][1][2];
+
+ c[2][2] = c4[2][2][2][2];
+ c[2][3] = c4[2][2][0][1];
+ c[2][4] = c4[2][2][0][2];
+ c[2][5] = c4[2][2][1][2];
+
+ c[3][3] = c4[0][1][0][1];
+ c[3][4] = c4[0][1][0][2];
+ c[3][5] = c4[0][1][1][2];
+
+ c[4][4] = c4[0][2][0][2];
+ c[4][5] = c4[0][2][1][2];
+
+ c[5][5] = c4[1][2][1][2];
+ /* fill in the lower half */
+ for(i=0;i<6;i++)
+ for(j=i+1;j<6;j++)
+ c[j][i] = c[i][j];
+}
+
+
+
+
+
+#endif
Modified: mc/3D/CitcomS/trunk/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Element_calculations.c 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/lib/Element_calculations.c 2010-09-13 07:13:07 UTC (rev 17185)
@@ -46,7 +46,7 @@
static void get_elt_tr(struct All_variables *, int , int , double [24], int );
static void get_elt_tr_pseudo_surf(struct All_variables *, int , int , double [24], int );
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
#include "anisotropic_viscosity.h"
#endif
@@ -302,7 +302,7 @@
const int ends = ENODES3D;
const int dims=E->mesh.nsd;
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
double D[VPOINTS3D+1][6][6],n[3],btmp[6];
int l1,l2;
#endif
@@ -318,11 +318,17 @@
for(k=1;k<=vpts;k++) {
off = (el-1)*vpts+k;
W[k]=g_point[k].weight[dims-1]*E->GDA[lev][m][el].vpt[k]*E->EVI[lev][m][off];
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
- if(E->viscosity.allow_orthotropic_viscosity){/* allow for a possibly anisotropic viscosity */
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ /* allow for a possibly anisotropic viscosity */
+
+ if(E->viscosity.allow_anisotropic_viscosity == 1){ /* orthotropic */
n[0] = E->EVIn1[lev][m][off];n[1] = E->EVIn2[lev][m][off];n[2] = E->EVIn3[lev][m][off]; /* Cartesian directors */
/* get the viscosity factor matrix and convert to CitcomS spherical */
get_constitutive_orthotropic_viscosity(D[k],E->EVI2[lev][m][off],n,TRUE,rtf[1][k],rtf[2][k]);
+ }else if(E->viscosity.allow_anisotropic_viscosity == 2){
+ /* transversely isotropic */
+ n[0] = E->EVIn1[lev][m][off];n[1] = E->EVIn2[lev][m][off];n[2] = E->EVIn3[lev][m][off]; /* Cartesian directors */
+ get_constitutive_ti_viscosity(D[k],E->EVI2[lev][m][off],0,n,TRUE,rtf[1][k],rtf[2][k]);
}
#endif
}
@@ -337,8 +343,8 @@
bdbmu[2][1]=bdbmu[2][2]=bdbmu[2][3]=
bdbmu[3][1]=bdbmu[3][2]=bdbmu[3][3]=0.0;
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
- if(E->viscosity.allow_orthotropic_viscosity){
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity){
for(i=1;i<=dims;i++)
for(j=1;j<=dims;j++)
for(k=1;k<=vpts;k++){ /* */
@@ -377,7 +383,7 @@
bdbmu[i][j] -= W[k] * two_thirds *
( ba[a][k][i][1] + ba[a][k][i][2] + ba[a][k][i][3] ) *
( ba[b][k][j][1] + ba[b][k][j][2] + ba[b][k][j][3] );
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
}
#endif
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-09-13 07:13:07 UTC (rev 17185)
@@ -44,9 +44,10 @@
#include "composition_related.h"
#include "element_definitions.h"
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
#include "anisotropic_viscosity.h"
#endif
+
#ifdef USE_GGRD
#include "hc.h" /* ggrd and hc packages */
@@ -1314,14 +1315,14 @@
}
-
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
/*
read in anisotropic viscosity from a directory which holds
-vis2.grd for the viscosity factors (1 - eta_S/eta)
+vis2.grd for the viscosity factors, read in log10(eta_S/eta)
nr.grd, nt.grd, np.grd for the directors
@@ -1334,7 +1335,7 @@
int m,el,i,j,k,l,inode,i1,i2,elxlz,elxlylz,ind,nel;
int llayer,nox,noy,noz,level,lselect,idim,elx,ely,elz;
char gmt_string[10],char_dummy;
- double vis2,ntheta,nphi,nr,rout[3],xloc[4],nlen;
+ double vis2,log_vis,ntheta,nphi,nr,rout[3],xloc[4],nlen;
double cvec[3],base[9];
char tfilename[1000];
static ggrd_boolean shift_to_pos_lon = FALSE;
@@ -1350,12 +1351,9 @@
elxlz = elx * elz;
elxlylz = elxlz * ely;
-#ifndef CITCOM_ALLOW_ORTHOTROPIC_VISC
- fprintf(stderr,"ggrd_read_anivisc_from_file: error, need to compile with CITCOM_ALLOW_ORTHOTROPIC_VISC\n");
- parallel_process_termination();
-#endif
- if(!E->viscosity.allow_orthotropic_viscosity)
- myerror(E,"ggrd_read_anivisc_from_file: called, but allow_orthotropic_viscosity is FALSE?!");
+
+ if(E->viscosity.allow_anisotropic_viscosity == 0)
+ myerror(E,"ggrd_read_anivisc_from_file: called, but allow_anisotropic_viscosity is FALSE?!");
/* isotropic default */
for(i=E->mesh.gridmin;i <= E->mesh.gridmax;i++){
@@ -1468,7 +1466,7 @@
/* vis2 */
if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],
- vis2_grd,&vis2,FALSE,shift_to_pos_lon)){
+ vis2_grd,&log_vis,FALSE,shift_to_pos_lon)){
fprintf(stderr,"ggrd_read_anivisc_from_file: interpolation error at lon: %g lat: %g\n",
rout[2]*180/M_PI,90-rout[1]*180/M_PI);
parallel_process_termination();
@@ -1509,6 +1507,8 @@
to
Cartesian */
convert_pvec_to_cvec_d(nr,ntheta,nphi,base,cvec);
+ /* transform */
+ vis2 = 1.0 - pow(10.0,log_vis);
for(l=1;l <= vpts;l++){ /* assign to all integration points */
ind = (el-1)*vpts + l;
E->EVI2[E->mesh.gridmax][m][ind] = vis2;
@@ -1532,6 +1532,7 @@
}
-#endif
+#endif /* for ANISOTROPIC */
+#endif /* for GGRD */
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2010-09-13 07:13:07 UTC (rev 17185)
@@ -1055,8 +1055,9 @@
} /* end for cap and i & j */
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
- if(E->viscosity.allow_orthotropic_viscosity){
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity){ /* any anisotropic
+ viscosity */
for(i=E->mesh.gridmin;i<=E->mesh.gridmax;i++)
for (j=1;j<=E->sphere.caps_per_proc;j++) {
nel = E->lmesh.NEL[i];
@@ -1078,9 +1079,11 @@
parallel_process_termination();
}
}
- E->viscosity.orthotropic_viscosity_init = FALSE;
+ E->viscosity.anisotropic_viscosity_init = FALSE;
if(E->parallel.me == 0)
- fprintf(stderr,"allocated for anisotropic viscosity levmax %i\n",E->mesh.gridmax);
+ fprintf(stderr,"allocated for anisotropic viscosity (%s) levmax %i\n",
+ (E->viscosity.allow_anisotropic_viscosity == 1)?("orthotropic"):("transversely isotropic"),
+ E->mesh.gridmax);
}
#endif
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2010-09-13 07:13:07 UTC (rev 17185)
@@ -58,7 +58,7 @@
extern void get_STD_topo(struct All_variables *, float**, float**,
float**, float**, int);
extern void get_CBF_topo(struct All_variables *, float**, float**);
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
#include "anisotropic_viscosity.h"
void output_avisc(struct All_variables *, int);
#endif
@@ -112,7 +112,7 @@
output_velo(E, cycles);
output_visc(E, cycles);
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
output_avisc(E, cycles);
#endif
@@ -322,14 +322,14 @@
return;
}
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
void output_avisc(struct All_variables *E, int cycles)
{
int i, j;
char output_file[255];
FILE *fp1;
int lev = E->mesh.levmax;
- if(E->viscosity.allow_orthotropic_viscosity){
+ if(E->viscosity.allow_anisotropic_viscosity){
sprintf(output_file,"%s.avisc.%d.%d", E->control.data_file,
E->parallel.me, cycles);
fp1 = output_open(output_file, "w");
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2010-09-13 07:13:07 UTC (rev 17185)
@@ -117,7 +117,7 @@
int open_file_zipped(char *, FILE **,struct All_variables *);
void gzip_file(char *);
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
#include "anisotropic_viscosity.h"
void gzdir_output_avisc(struct All_variables *, int);
#endif
@@ -167,7 +167,7 @@
else new VTK output won't
work */
gzdir_output_visc(E, out_cycles);
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
gzdir_output_avisc(E, out_cycles);
#endif
@@ -760,6 +760,8 @@
return;
}
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+
/*
anisotropic viscosity
*/
@@ -775,7 +777,7 @@
MPI_Status mpi_stat;
int mpi_rc;
int mpi_inmsg, mpi_success_message = 1;
- if(E->viscosity.allow_orthotropic_viscosity){
+ if(E->viscosity.allow_anisotropic_viscosity){
if(E->output.gzdir.vtk_io < 2){
snprintf(output_file,255,
@@ -820,8 +822,8 @@
return;
}
+#endif
-
void gzdir_output_surf_botm(struct All_variables *E, int cycles)
{
int i, j, s;
Modified: mc/3D/CitcomS/trunk/lib/Solver_multigrid.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Solver_multigrid.c 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/lib/Solver_multigrid.c 2010-09-13 07:13:07 UTC (rev 17185)
@@ -253,8 +253,8 @@
viscD[m]=(float *)malloc((1+vpts*E->lmesh.NEL[lv-1])*sizeof(float));
}
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC /* allow for anisotropy */
- if(E->viscosity.allow_orthotropic_viscosity){
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC /* allow for anisotropy */
+ if(E->viscosity.allow_anisotropic_viscosity){
for(lv=E->mesh.levmax;lv>E->mesh.levmin;lv--) {
sl_minus = lv -1;
if (E->viscosity.smooth_cycles==1) {
@@ -339,7 +339,7 @@
}
*/
}
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
}
#endif
for(m=1;m<=E->sphere.caps_per_proc;m++) {
Modified: mc/3D/CitcomS/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2010-09-13 07:13:07 UTC (rev 17185)
@@ -48,7 +48,7 @@
float** , float** , float** ,
float** , float** );
void stress_conform_bcs(struct All_variables *);
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
#include "anisotropic_viscosity.h"
#endif
@@ -229,7 +229,7 @@
double pre[9],tww[9],rtf[4][9];
double velo_scaling, stress_scaling, mass_fac;
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
double D[6][6],n[3],eps[6],str[6];
#endif
struct Shape_function_dA *dOmega;
@@ -365,8 +365,8 @@
for(j=1;j<=vpts;j++)
dilation[j] = (Vxyz[1][j] + Vxyz[2][j] + Vxyz[3][j]) / 3.0;
}
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
- if(E->viscosity.allow_orthotropic_viscosity){
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity){ /* general anisotropic */
for(j=1;j <= vpts;j++) {
l1 = (e-1)*vpts+j;
n[0] = E->EVIn1[E->mesh.levmax][m][l1]; /* Cartesian directors */
@@ -377,7 +377,10 @@
CitcomS convection
*/
- get_constitutive_orthotropic_viscosity(D,E->EVI2[E->mesh.levmax][m][l1],n,TRUE,rtf[1][j],rtf[2][j]);
+ if(E->viscosity.allow_anisotropic_viscosity == 1)
+ get_constitutive_orthotropic_viscosity(D,E->EVI2[E->mesh.levmax][m][l1],n,TRUE,rtf[1][j],rtf[2][j]);
+ else if(E->viscosity.allow_anisotropic_viscosity == 2)
+ get_constitutive_ti_viscosity(D,E->EVI2[E->mesh.levmax][m][l1],0,n,TRUE,rtf[1][j],rtf[2][j]);
/* deviatoric stress, pressure will be added later */
eps[0] = Vxyz[1][j] - dilation[j]; /* strain-rates */
@@ -413,7 +416,7 @@
div += Vxyz[7][j]*dOmega->vpt[j]; /* divergence */
vor += Vxyz[8][j]*dOmega->vpt[j]; /* vorticity */
}
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
}
#endif
/* normalize by volume */
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2010-09-13 07:13:07 UTC (rev 17185)
@@ -48,7 +48,7 @@
static void low_viscosity_channel_factor(struct All_variables *E, float *F);
static void low_viscosity_wedge_factor(struct All_variables *E, float *F);
void parallel_process_termination();
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
#include "anisotropic_viscosity.h"
#endif
@@ -65,18 +65,18 @@
input_string("visc_layer_file", E->viscosity.layer_file,"visc.dat",m);
- input_boolean("allow_orthotropic_viscosity",&(E->viscosity.allow_orthotropic_viscosity),"off",m);
-#ifndef CITCOM_ALLOW_ORTHOTROPIC_VISC
- if(E->viscosity.allow_orthotropic_viscosity){ /* error */
- fprintf(stderr,"error: allow_orthotropic_viscosity is set, but code not compiled with CITCOM_ALLOW_ORTHOTROPIC_VISC\n");
+ input_int("allow_anisotropic_viscosity",&(E->viscosity.allow_anisotropic_viscosity),"0",m);
+#ifndef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity){ /* error */
+ fprintf(stderr,"error: allow_anisotropic_viscosity is not zero, but code not compiled with CITCOM_ALLOW_ANISOTROPIC_VISC\n");
parallel_process_termination();
}
#else
- if(E->viscosity.allow_orthotropic_viscosity){ /* read additional
+ if(E->viscosity.allow_anisotropic_viscosity){ /* read additional
parameters for
anisotropic
viscosity */
- input_int("anisotropic_init",&(E->viscosity.anisotropic_init),"0",m);
+ input_int("anisotropic_init",&(E->viscosity.anisotropic_init),"0",m); /* 0: isotropic 1: random 2: read in director and log10(eta_s/eta) */
input_string("anisotropic_init_dir",(E->viscosity.anisotropic_init_dir),"",m); /* directory
for
ggrd
@@ -275,9 +275,9 @@
double *TG;
const int vpts = vpoints[E->mesh.nsd];
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
- if(E->viscosity.allow_orthotropic_viscosity){
- if(!E->viscosity.orthotropic_viscosity_init)
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity){
+ if(!E->viscosity.anisotropic_viscosity_init)
set_anisotropic_viscosity_at_element_level(E,1);
else
set_anisotropic_viscosity_at_element_level(E,0);
@@ -345,8 +345,8 @@
visc_from_nodes_to_gint(E,visc,evisc,E->mesh.levmax);
}
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC /* allow for anisotropy */
- if(E->viscosity.allow_orthotropic_viscosity){
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC /* allow for anisotropy */
+ if(E->viscosity.allow_anisotropic_viscosity){
visc_from_gint_to_nodes(E,E->EVI2[E->mesh.levmax], E->VI2[E->mesh.levmax],E->mesh.levmax);
visc_from_gint_to_nodes(E,E->EVIn1[E->mesh.levmax], E->VIn1[E->mesh.levmax],E->mesh.levmax);
visc_from_gint_to_nodes(E,E->EVIn2[E->mesh.levmax], E->VIn2[E->mesh.levmax],E->mesh.levmax);
@@ -354,6 +354,7 @@
normalize_director_at_nodes(E,E->VIn1[E->mesh.levmax],E->VIn2[E->mesh.levmax],E->VIn3[E->mesh.levmax],E->mesh.levmax);
if(E->viscosity.SMOOTH){
+ if(E->parallel.me == 0)fprintf(stderr,"WARNING: smoothing anisotropic viscosity, perhaps not a good idea\n");
visc_from_nodes_to_gint(E,E->VI2[E->mesh.levmax],E->EVI2[E->mesh.levmax],E->mesh.levmax);
visc_from_nodes_to_gint(E,E->VIn1[E->mesh.levmax],E->EVIn1[E->mesh.levmax],E->mesh.levmax);
visc_from_nodes_to_gint(E,E->VIn2[E->mesh.levmax],E->EVIn2[E->mesh.levmax],E->mesh.levmax);
Modified: mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h 2010-09-13 07:13:07 UTC (rev 17185)
@@ -1,20 +1,28 @@
#ifndef __CITCOM_READ_ANIVISC_HEADER__
+void get_constitutive_ti_viscosity(double [6][6], double, double ,double [3], int ,
+ double , double) ;
-void get_constitutive_orthotropic_viscosity(double [6][6], double,double [3], int, double, double) ;
+void get_constitutive_orthotropic_viscosity(double [6][6], double,
+ double [3], int,
+ double , double ) ;
+void set_anisotropic_viscosity_at_element_level(struct All_variables *, int ) ;
-void set_anisotropic_viscosity_at_element_level(struct All_variables *,int);
-void normalize_director_at_nodes(struct All_variables *,float **,float **, float **, int );
-void normalize_director_at_gint(struct All_variables *,float **,float **, float **, int );
+void normalize_director_at_nodes(struct All_variables *, float **, float **, float **, int);
+void normalize_director_at_gint(struct All_variables *, float **, float **, float **, int);
+void conv_cart4x4_to_spherical(double [3][3][3][3], double, double, double [3][3][3][3]);
+void conv_cart6x6_to_spherical(double [6][6], double, double, double [6][6]);
+void get_citcom_spherical_rot(double, double, double [3][3]);
+void get_orth_delta(double [3][3][3][3], double [3]);
+void rot_4x4(double [3][3][3][3], double [3][3], double [3][3][3][3]);
+void copy_4x4(double [3][3][3][3], double [3][3][3][3]);
+void copy_6x6(double [6][6], double [6][6]);
+void c4fromc6(double [3][3][3][3], double [6][6]);
+void c6fromc4(double [6][6], double [3][3][3][3]);
+void print_6x6_mat(FILE *, double [6][6]);
+void zero_6x6(double [6][6]);
+void zero_4x4(double [3][3][3][3]);
+void rotate_ti6x6_to_director(double [6][6],double [3]);
-
-
-
-
-void conv_cart4x4_to_spherical(double [3][3][3][3],
- double , double, double [3][3][3][3]);
-void get_delta(double [3][3][3][3],double [3]);
-void rot_4x4(double [3][3][3][3],double [3][3], double [3][3][3][3]);
-void print_6x6_mat(FILE *, double [6][6]);
#define __CITCOM_READ_ANIVISC_HEADER__
#endif
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2010-09-13 07:13:07 UTC (rev 17185)
@@ -784,7 +784,7 @@
float *Vi[NCS],*EVi[NCS];
float *VI[MAX_LEVELS][NCS],*EVI[MAX_LEVELS][NCS];
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
float *VI2[MAX_LEVELS][NCS],*EVI2[MAX_LEVELS][NCS];
float *VIn1[MAX_LEVELS][NCS],*EVIn1[MAX_LEVELS][NCS];
float *VIn2[MAX_LEVELS][NCS],*EVIn2[MAX_LEVELS][NCS];
Modified: mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2010-09-11 22:12:43 UTC (rev 17184)
+++ mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2010-09-13 07:13:07 UTC (rev 17185)
@@ -35,8 +35,8 @@
int SMOOTH;
int smooth_cycles;
- int allow_orthotropic_viscosity,orthotropic_viscosity_init;
-#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ int allow_anisotropic_viscosity,anisotropic_viscosity_init;
+#ifdef CITCOM_ALLOW_ANISOTROPIC_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 */
More information about the CIG-COMMITS
mailing list