[cig-commits] r17187 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Mon Sep 13 07:07:10 PDT 2010
Author: becker
Date: 2010-09-13 07:07:10 -0700 (Mon, 13 Sep 2010)
New Revision: 17187
Modified:
mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
mc/3D/CitcomS/trunk/lib/Drive_solvers.c
mc/3D/CitcomS/trunk/lib/Element_calculations.c
mc/3D/CitcomS/trunk/lib/Ggrd_handling.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:
Added option to start with isotropic viscosity during first anisotropic
iteration.
Modified: mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c 2010-09-13 12:01:43 UTC (rev 17186)
+++ mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c 2010-09-13 14:07:10 UTC (rev 17187)
@@ -51,6 +51,37 @@
#define CITCOM_DELTA(i,j) ((i==j)?(1.0):(0.0))
+void get_constitutive(double D[6][6], int lev, int m,
+ int off, double theta, double phi,
+ struct All_variables *E)
+{
+ double n[3];
+ const int convert_to_spherical = TRUE;
+ if(E->viscosity.allow_anisotropic_viscosity){
+ if((E->monitor.solution_cycles == 0)&&
+ (E->viscosity.anivisc_start_from_iso)&&(E->monitor.visc_iter_count == 0))
+ /* first iteration of "stress-dependent" loop for first timestep */
+ get_constitutive_isotropic(D);
+ else{
+ /* allow for a possibly anisotropic viscosity */
+ n[0] = E->EVIn1[lev][m][off];
+ n[1] = E->EVIn2[lev][m][off];
+ n[2] = E->EVIn3[lev][m][off]; /* Cartesian directors */
+ if(E->viscosity.allow_anisotropic_viscosity == 1){ /* orthotropic */
+ get_constitutive_orthotropic_viscosity(D,E->EVI2[lev][m][off],
+ n,convert_to_spherical,theta,phi);
+ }else if(E->viscosity.allow_anisotropic_viscosity == 2){
+ /* transversely isotropic */
+ get_constitutive_ti_viscosity(D,E->EVI2[lev][m][off],0.,
+ n,convert_to_spherical,theta,phi);
+ }
+ }
+ }else{
+ get_constitutive_isotropic(D);
+ }
+}
+
+
/*
transversely isotropic viscosity following Han and Wahr
@@ -86,14 +117,7 @@
double nlen,delta_vis2;
int ani;
/* 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 */
-
+ get_constitutive_isotropic(D);
ani = FALSE;
if((fabs(delta_vis) > 3e-15) || (fabs(gamma_vis) > 3e-15)){
ani = TRUE;
@@ -148,15 +172,8 @@
double delta[3][3][3][3];
int ani;
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 */
- D[3][3] = 1.0; /* xy = tp */
- D[4][4] = 1.0; /* xz = rt */
- D[5][5] = 1.0; /* yz = rp */
-
+ /* start with isotropic */
+ get_constitutive_isotropic(D);
/* get Cartesian anisotropy matrix */
if(fabs(delta_vis) > 3e-15){
ani = TRUE;
@@ -217,7 +234,17 @@
//print_6x6_mat(stderr,D);fprintf(stderr,"\n");
}
}
-
+void get_constitutive_isotropic(double D[6][6])
+{
+ /* 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 */
+}
void set_anisotropic_viscosity_at_element_level(struct All_variables *E, int init_stage)
{
int i,j,k,l,off,nel;
@@ -287,7 +314,7 @@
fprintf(stderr,"set_anisotropic_viscosity_at_element_level: anisotropic_init mode 2 requires USE_GGRD compilation\n");
parallel_process_termination();
#endif
- ggrd_read_anivisc_from_file(E,1);
+ ggrd_read_anivisc_from_file(E,(E->sphere.caps == 12)?(TRUE):(FALSE));
break;
default:
fprintf(stderr,"set_anisotropic_viscosity_at_element_level: anisotropic_init %i undefined\n",
@@ -352,7 +379,8 @@
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])
+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);
@@ -370,27 +398,31 @@
*/
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],
+ double a[3][3][3][3],b[3][3][3][3],rot[3][3],test[3],testr[3],prod[3],
hlen2,x2,y2,xy,zm1;
/* normalize */
normalize_vec3d((n+0),(n+1),(n+2));
/* calc aux variable */
x2 = n[0]*n[0];y2 = n[1]*n[1];xy = n[0]*n[1];
hlen2 = x2 + y2;zm1 = n[2]-1;
- /* rotation matrix to get {0,0,1} to {x,y,z} */
- rot[0][0] = (y2 + x2*n[2])/hlen2;
- rot[0][1] = (xy*zm1)/hlen2;
- rot[0][2] = n[0];
- rot[1][0] = rot[0][1];
- rot[1][1] = (x2 + y2*n[2])/hlen2;
- rot[1][2] = n[1];
- rot[2][0] = -n[0];
- rot[2][1] = -n[1];
- rot[2][2] = n[2];
- /* rotate the D matrix */
- c4fromc6(a,D);
- rot_4x4(a,rot,b);
- c6fromc4(D,b);
+ if(hlen2 > 3e-15){
+ /* rotation matrix to get {0,0,1} to {x,y,z} */
+ rot[0][0] = (y2 + x2*n[2])/hlen2;
+ rot[0][1] = (xy*zm1)/hlen2;
+ rot[0][2] = n[0];
+ rot[1][0] = rot[0][1];
+ rot[1][1] = (x2 + y2*n[2])/hlen2;
+ rot[1][2] = n[1];
+ rot[2][0] = -n[0];
+ rot[2][1] = -n[1];
+ rot[2][2] = n[2];
+
+ /* rotate the D matrix */
+ c4fromc6(a,D);
+ rot_4x4(a,rot,b);
+ c6fromc4(D,b);
+ } /* already oriented right */
+
}
void get_citcom_spherical_rot(double theta, double phi, double rot[3][3]){
@@ -597,7 +629,15 @@
c[j][i] = c[i][j];
}
-
+void mat_mult_vec_3x3(double a[3][3],double b[3],double c[3])
+{
+ int i,j;
+ for(i=0;i<3;i++){
+ c[i]=0;
+ for(j=0;j<3;j++)
+ c[i] += a[i][j] * b[j];
+ }
+}
void normalize_vec3(float *x, float *y, float *z)
{
double len = 0.;
@@ -614,6 +654,11 @@
len = sqrt(len);
*x /= len;*y /= len;*z /= len;
}
+void cross_product(double a[3],double b[3],double c[3])
+{
+ c[0]=a[1]*b[2]-a[2]*b[1];
+ c[1]=a[2]*b[0]-a[0]*b[2];
+ c[2]=a[0]*b[1]-a[1]*b[0];
+}
-
#endif
Modified: mc/3D/CitcomS/trunk/lib/Drive_solvers.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2010-09-13 12:01:43 UTC (rev 17186)
+++ mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2010-09-13 14:07:10 UTC (rev 17187)
@@ -34,9 +34,9 @@
double global_vdot();
double vnorm_nonnewt();
int need_visc_update(struct All_variables *);
+int need_to_iterate(struct All_variables *);
-
/************************************************************/
void general_stokes_solver_setup(struct All_variables *E)
@@ -68,12 +68,14 @@
void remove_rigid_rot();
double Udot_mag, dUdot_mag;
- int m,count,i;
+ int m,i;
double *oldU[NCS], *delta_U[NCS];
const int neq = E->lmesh.neq;
+ E->monitor.visc_iter_count = 0; /* first solution */
+
velocities_conform_bcs(E,E->U);
assemble_forces(E,0);
@@ -81,10 +83,11 @@
get_system_viscosity(E,1,E->EVI[E->mesh.levmax],E->VI[E->mesh.levmax]);
construct_stiffness_B_matrix(E);
}
-
+
solve_constrained_flow_iterative(E);
- if (E->viscosity.SDEPV || E->viscosity.PDEPV) {
+
+ if (need_to_iterate(E)) {
/* outer iterations for velocity dependent viscosity */
for (m=1;m<=E->sphere.caps_per_proc;m++) {
@@ -95,8 +98,8 @@
}
Udot_mag=dUdot_mag=0.0;
- count=1;
+ E->monitor.visc_iter_count++;
while (1) {
@@ -112,19 +115,20 @@
if(E->parallel.me==0){
fprintf(stderr,"Stress dep. visc./plast.: DUdot = %.4e (%.4e) for iteration %d\n",
- dUdot_mag,Udot_mag,count);
+ dUdot_mag,Udot_mag,E->monitor.visc_iter_count);
fprintf(E->fp,"Stress dep. visc./plast.: DUdot = %.4e (%.4e) for iteration %d\n",
- dUdot_mag,Udot_mag,count);
+ dUdot_mag,Udot_mag,E->monitor.visc_iter_count);
fflush(E->fp);
}
- if ((count>50) || (dUdot_mag < E->viscosity.sdepv_misfit))
+ if ((E->monitor.visc_iter_count > 50) ||
+ (dUdot_mag < E->viscosity.sdepv_misfit))
break;
get_system_viscosity(E,1,E->EVI[E->mesh.levmax],E->VI[E->mesh.levmax]);
construct_stiffness_B_matrix(E);
solve_constrained_flow_iterative(E);
- count++;
+ E->monitor.visc_iter_count++;
} /*end while*/
@@ -133,7 +137,7 @@
free((void *) delta_U[m]);
}
- } /*end if SDEPV or PDEPV */
+ } /*end if we need iterations */
/* remove the rigid rotation component from the velocity solution */
if((E->sphere.caps == 12) &&
@@ -167,7 +171,23 @@
}
}
}
-
+int need_to_iterate(struct All_variables *E){
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ /* anisotropic viscosity */
+ if(E->viscosity.allow_anisotropic_viscosity){
+ if((E->monitor.solution_cycles == 0) &&
+ E->viscosity.anivisc_start_from_iso) /* first step will be
+ solved isotropically at
+ first */
+ return TRUE;
+ else
+ return (E->viscosity.SDEPV || E->viscosity.PDEPV)?(TRUE):(FALSE);
+ }
+#else
+ /* regular operation */
+ return ((E->viscosity.SDEPV || E->viscosity.PDEPV)?(TRUE):(FALSE));
+#endif
+}
void general_stokes_solver_pseudo_surf(struct All_variables *E)
{
void solve_constrained_flow_iterative();
Modified: mc/3D/CitcomS/trunk/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Element_calculations.c 2010-09-13 12:01:43 UTC (rev 17186)
+++ mc/3D/CitcomS/trunk/lib/Element_calculations.c 2010-09-13 14:07:10 UTC (rev 17187)
@@ -303,7 +303,7 @@
const int dims=E->mesh.nsd;
#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
- double D[VPOINTS3D+1][6][6],n[3],btmp[6];
+ double D[VPOINTS3D+1][6][6],btmp[6];
int l1,l2;
#endif
@@ -320,16 +320,7 @@
W[k]=g_point[k].weight[dims-1]*E->GDA[lev][m][el].vpt[k]*E->EVI[lev][m][off];
#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]);
- }
+ get_constitutive(D[k],lev,m,off,rtf[1][k],rtf[2][k],E);
#endif
}
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-09-13 12:01:43 UTC (rev 17186)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-09-13 14:07:10 UTC (rev 17187)
@@ -1368,11 +1368,6 @@
}
}
}
- /*
-
- rest
-
- */
if(is_global) /* decide on GMT flag */
sprintf(gmt_string,GGRD_GMT_GLOBAL_STRING); /* global */
else
Modified: mc/3D/CitcomS/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2010-09-13 12:01:43 UTC (rev 17186)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2010-09-13 14:07:10 UTC (rev 17187)
@@ -369,18 +369,12 @@
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 */
- n[1] = E->EVIn2[E->mesh.levmax][m][l1];
- n[2] = E->EVIn3[E->mesh.levmax][m][l1];
/*
get viscosity matrix and convert to spherical system in
CitcomS convection
*/
- 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]);
+ get_constitutive(D,E->mesh.levmax,m,l1,rtf[1][j],rtf[2][j],E);
/* deviatoric stress, pressure will be added later */
eps[0] = Vxyz[1][j] - dilation[j]; /* strain-rates */
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2010-09-13 12:01:43 UTC (rev 17186)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2010-09-13 14:07:10 UTC (rev 17187)
@@ -86,6 +86,11 @@
<0: assign to layer = anivisc_layer
*/
+ input_boolean("anivisc_start_from_iso",
+ &(E->viscosity.anivisc_start_from_iso),"on",m); /* start
+ from
+ isotropic
+ solution? */
}
#endif
Modified: mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h 2010-09-13 12:01:43 UTC (rev 17186)
+++ mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h 2010-09-13 14:07:10 UTC (rev 17187)
@@ -24,6 +24,11 @@
void rotate_ti6x6_to_director(double [6][6],double [3]);
void normalize_vec3(float *, float *, float *);
void normalize_vec3d(double *, double *, double *);
+void mat_mult_vec_3x3(double [3][3],double [3],double [3]);
+void cross_product(double [3],double [3],double [3]);
+void get_constitutive_isotropic(double [6][6]);
+void get_constitutive(double [6][6], int , int , int , double , double ,
+ struct All_variables *);
#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-13 12:01:43 UTC (rev 17186)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2010-09-13 14:07:10 UTC (rev 17187)
@@ -378,6 +378,8 @@
int solution_cycles;
int solution_cycles_init;
+ int visc_iter_count;
+
int stop_topo_loop;
int topo_loop;
Modified: mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2010-09-13 12:01:43 UTC (rev 17186)
+++ mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2010-09-13 14:07:10 UTC (rev 17187)
@@ -37,8 +37,9 @@
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_start_from_iso; /* start from isotropic solution? */
+ 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
More information about the CIG-COMMITS
mailing list