[cig-commits] r17197 - mc/3D/CitcomCU/trunk/src
becker at geodynamics.org
becker at geodynamics.org
Wed Sep 15 22:57:34 PDT 2010
Author: becker
Date: 2010-09-15 22:57:34 -0700 (Wed, 15 Sep 2010)
New Revision: 17197
Added:
mc/3D/CitcomCU/trunk/src/Anisotropic_viscosity.c
mc/3D/CitcomCU/trunk/src/Makefile.gzdir_ani
mc/3D/CitcomCU/trunk/src/anisotropic_viscosity.h
Modified:
mc/3D/CitcomCU/trunk/src/Citcom.c
mc/3D/CitcomCU/trunk/src/Construct_arrays.c
mc/3D/CitcomCU/trunk/src/Convection.c
mc/3D/CitcomCU/trunk/src/Drive_solvers.c
mc/3D/CitcomCU/trunk/src/Element_calculations.c
mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
mc/3D/CitcomCU/trunk/src/Instructions.c
mc/3D/CitcomCU/trunk/src/Makefile.gzdir
mc/3D/CitcomCU/trunk/src/Output_gzdir.c
mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c
mc/3D/CitcomCU/trunk/src/Parallel_related.c
mc/3D/CitcomCU/trunk/src/Solver_multigrid.c
mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c
mc/3D/CitcomCU/trunk/src/Topo_gravity.c
mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
mc/3D/CitcomCU/trunk/src/global_defs.h
mc/3D/CitcomCU/trunk/src/prototypes.h
mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h
Log:
Added facility to read in more than four material layers by means of
z_layer and r_layer keyword for cart and spherical, resp. default uses
four layers and the z_... and r_... settings where ... is lith, 410,
and lmantle (i.e. backward compatible), same as in CitcomS
Added consistent initialization of the random number generator.
First time step computation time is now computed right (no division by zero).
The Uzawa iteration log message has been modified to print absolute
velocity and pressure values, the latter are now also stored.
Test implementation of anisotropic viscosity, invisible to regular
compile.
Added: mc/3D/CitcomCU/trunk/src/Anisotropic_viscosity.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Anisotropic_viscosity.c (rev 0)
+++ mc/3D/CitcomCU/trunk/src/Anisotropic_viscosity.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -0,0 +1,625 @@
+/*
+
+ orthotropic viscosity following Muehlhaus, Moresi, Hobbs and Dufour
+ (PAGEOPH, 159, 2311, 2002)
+
+ tranverse isotropy following Han and Wahr (PEPI, 102, 33, 1997)
+
+
+*/
+
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+
+#include <math.h>
+#include "element_definitions.h"
+#include "global_defs.h"
+#include "anisotropic_viscosity.h"
+void calc_cbase_at_tp(float , float , float *);
+void calc_cbase_at_tp_d(double , double , double *);
+
+#define CITCOM_DELTA(i,j) ((i==j)?(1.0):(0.0))
+void ggrd_read_anivisc_from_file(struct All_variables *);
+
+void get_constitutive(double D[6][6], int lev,
+ int off, double theta, double phi,
+ int convert_to_spherical,
+ struct All_variables *E)
+{
+ double n[3];
+ 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][off];
+ n[1] = E->EVIn2[lev][off];
+ n[2] = E->EVIn3[lev][off]; /* Cartesian directors */
+ if(E->viscosity.allow_anisotropic_viscosity == 1){ /* orthotropic */
+ get_constitutive_orthotropic_viscosity(D,E->EVI2[lev][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][off],0.,n,convert_to_spherical,theta,phi);
+ }
+ }
+ }else{
+ get_constitutive_isotropic(D);
+ }
+}
+
+
+/*
+
+transversely isotropic viscosity following Han and Wahr
+
+
+\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 direction) and will
+normalize n, if not already normalized
+
+
+*/
+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;
+ /* isotropic part, in units of iso_visc */
+ get_constitutive_isotropic(D);
+ ani = FALSE;
+ if((fabs(delta_vis) > 3e-15) || (fabs(gamma_vis) > 3e-15)){
+ ani = TRUE;
+ /* get Cartesian anisotropy matrix by adding anisotropic
+ components */
+ D[0][0] += gamma_vis;
+ D[1][0] = D[0][1] = gamma_vis;
+ D[1][1] = D[0][0];
+ D[2][2] += 2.*gamma_vis;
+ D[4][4] -= delta_vis;
+ D[5][5] = D[4][4];
+ /*
+ the rotation routine takes care of normalization and will normalize n
+ */
+ //print_6x6_mat(stderr,D);
+ rotate_ti6x6_to_director(D,n); /* rotate such that the generic z
+ preferred axis is aligned with
+ the director */
+ //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, specify in cartesian
+
+
+ where delta_vis = (1 - eta_S/eta)
+
+
+ */
+void get_constitutive_orthotropic_viscosity(double D[6][6], double delta_vis,
+ double n[3], int convert_to_spherical,
+ double theta, double phi)
+{
+ double nlen,delta_vis2;
+ double delta[3][3][3][3];
+ int ani;
+ ani=FALSE;
+ /* start with isotropic */
+ get_constitutive_isotropic(D);
+ /* 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)
+
+ this routine will normalize n, just in case
+
+ */
+ delta_vis2 = 2.0*delta_vis;
+ /* s_xx = s_tt */
+ D[0][0] -= delta_vis2 * delta[0][0][0][0]; /* * e_xx */
+ D[0][1] -= delta_vis2 * delta[0][0][1][1];
+ D[0][2] -= delta_vis2 * delta[0][0][2][2];
+ D[0][3] -= delta_vis * (delta[0][0][0][1]+delta[0][0][1][0]);
+ D[0][4] -= delta_vis * (delta[0][0][0][2]+delta[0][0][2][0]);
+ D[0][5] -= delta_vis * (delta[0][0][1][2]+delta[0][0][2][1]);
+
+ D[1][0] -= delta_vis2 * delta[1][1][0][0]; /* s_yy = s_pp */
+ D[1][1] -= delta_vis2 * delta[1][1][1][1];
+ D[1][2] -= delta_vis2 * delta[1][1][2][2];
+ D[1][3] -= delta_vis * (delta[1][1][0][1]+delta[1][1][1][0]);
+ D[1][4] -= delta_vis * (delta[1][1][0][2]+delta[1][1][2][0]);
+ D[1][5] -= delta_vis * (delta[1][1][1][2]+delta[1][1][2][1]);
+
+ D[2][0] -= delta_vis2 * delta[2][2][0][0]; /* s_zz = s_rr */
+ D[2][1] -= delta_vis2 * delta[2][2][1][1];
+ D[2][2] -= delta_vis2 * delta[2][2][2][2];
+ D[2][3] -= delta_vis * (delta[2][2][0][1]+delta[2][2][1][0]);
+ D[2][4] -= delta_vis * (delta[2][2][0][2]+delta[2][2][2][0]);
+ D[2][5] -= delta_vis * (delta[2][2][1][2]+delta[2][2][2][1]);
+
+ D[3][0] -= delta_vis2 * delta[0][1][0][0]; /* s_xy = s_tp */
+ D[3][1] -= delta_vis2 * delta[0][1][1][1];
+ D[3][2] -= delta_vis2 * delta[0][1][2][2];
+ D[3][3] -= delta_vis * (delta[0][1][0][1]+delta[0][1][1][0]);
+ D[3][4] -= delta_vis * (delta[0][1][0][2]+delta[0][1][2][0]);
+ D[3][5] -= delta_vis * (delta[0][1][1][2]+delta[0][1][2][1]);
+
+ D[4][0] -= delta_vis2 * delta[0][2][0][0]; /* s_xz = s_tr */
+ D[4][1] -= delta_vis2 * delta[0][2][1][1];
+ D[4][2] -= delta_vis2 * delta[0][2][2][2];
+ D[4][3] -= delta_vis * (delta[0][2][0][1]+delta[0][2][1][0]);
+ D[4][4] -= delta_vis * (delta[0][2][0][2]+delta[0][2][2][0]);
+ D[4][5] -= delta_vis * (delta[0][2][1][2]+delta[0][2][2][1]);
+
+ D[5][0] -= delta_vis2 * delta[1][2][0][0]; /* s_yz = s_pr */
+ D[5][1] -= delta_vis2 * delta[1][2][1][1];
+ D[5][2] -= delta_vis2 * delta[1][2][2][2];
+ D[5][3] -= delta_vis * (delta[1][2][0][1]+delta[1][2][1][0]);
+ 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 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,k,l,off,nel;
+ double vis2,n[3],u,v,s,r;
+ const int vpts = vpoints[E->mesh.nsd];
+
+ if(E->viscosity.allow_anisotropic_viscosity){
+ if(init_stage){
+ if(E->viscosity.anisotropic_viscosity_init)
+ myerror("anisotropic viscosity should not be initialized twice?!",E);
+ /* first call */
+ /* initialize anisotropic viscosity at element level, nodes will
+ get assigned later */
+ switch(E->viscosity.anisotropic_init){
+ case 0: /* isotropic */
+ if(E->parallel.me == 0)fprintf(stderr,"set_anisotropic_viscosity_at_element_level: initializing isotropic viscosity\n");
+ for(i=E->mesh.levmin;i <= E->mesh.levmax;i++){
+ nel = E->lmesh.NEL[i];
+ for(k=1;k <= nel;k++){
+ for(l=1;l <= vpts;l++){ /* assign to all integration points */
+ off = (k-1)*vpts + l;
+ E->EVI2[i][off] = 0.0;
+ E->EVIn1[i][off] = 1.0; E->EVIn2[i][off] = E->EVIn3[i][off] = 0.0;
+ }
+ }
+ }
+ break;
+ case 1: /*
+ random fluctuations, for testing a
+ worst case scenario
+
+ */
+ if(E->parallel.me == 0)fprintf(stderr,"set_anisotropic_viscosity_at_element_level: initializing random viscosity\n");
+ for(i=E->mesh.levmin;i <= E->mesh.levmax;i++){
+ nel = E->lmesh.NEL[i];
+ for(k=1;k <= nel;k++){
+ /* by element (srand48 call should be in citcom.c or somewhere? */
+ vis2 = drand48()*0.9; /* random fluctuation,
+ corresponding to same strength
+ (0) and 10 fold reduction
+ (0.9) */
+ /* get random vector */
+ do{
+ u = -1 + drand48()*2;v = -1 + drand48()*2;
+ s = u*u + v*v;
+ }while(s > 1);
+ r = 2.0 * sqrt(1.0-s );
+ n[0] = u * r; /* x */
+ n[1] = v * r; /* y */
+ n[2] = 2.0*s -1 ; /* z */
+ for(l=1;l <= vpts;l++){ /* assign to all integration points */
+ off = (k-1)*vpts + l;
+ E->EVI2[i][off] = vis2;
+ E->EVIn1[i][off] = n[0];
+ E->EVIn2[i][off] = n[1];
+ E->EVIn3[i][off] = n[2];
+ }
+ }
+ }
+ break;
+ case 2: /* from file */
+#ifndef USE_GGRD
+ 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);
+ break;
+ default:
+ fprintf(stderr,"set_anisotropic_viscosity_at_element_level: anisotropic_init %i undefined\n",
+ E->viscosity.anisotropic_init);
+ parallel_process_termination();
+ break;
+ }
+ E->viscosity.anisotropic_viscosity_init = TRUE;
+ /* end initialization stage */
+ }else{
+ /* standard operation every time step */
+
+
+ }
+ } /* end anisotropic viscosity branch */
+}
+
+
+void normalize_director_at_nodes(struct All_variables *E,float *n1,float *n2, float *n3, int lev)
+{
+ int n;
+ for(n=1;n<=E->lmesh.NNO[lev];n++){
+ normalize_vec3(&(n1[n]),&(n2[n]),&(n3[n]));
+ }
+}
+void normalize_director_at_gint(struct All_variables *E,float *n1,float *n2, float *n3, int lev)
+{
+ int e,i,off;
+ const int nsd=E->mesh.nsd;
+ const int vpts=vpoints[nsd];
+ for(e=1;e<=E->lmesh.NEL[lev];e++)
+ for(i=1;i<=vpts;i++) {
+ off = (e-1)*vpts+i;
+ normalize_vec3(&(n1[off]),&(n2[off]),&(n3[off]));
+ }
+}
+/*
+
+convert cartesian fourth order tensor (input c) to spherical, CitcomS
+format (output p)
+
+c and p cannot be the same matrix
+
+1: t 2: p 3: r
+
+(E only passed for debugging)
+
+*/
+
+ 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];
+ 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
+
+n will be normalized, just in case
+
+*/
+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],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;
+ 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]){
+ float base[9];
+ calc_cbase_at_tp((float)theta,(float)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]);
+}
+/*
+
+
+get fourth order anisotropy tensor for orthotropic viscosity from
+Muehlhaus et al. (2002)
+
+*/
+void get_orth_delta(double d[3][3][3][3],double n[3])
+{
+ int i,j,k,l;
+ double tmp;
+ normalize_vec3d((n+0),(n+1),(n+2));
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++)
+ for(k=0;k<3;k++)
+ for(l=0;l<3;l++){ /* eq. (4) from Muehlhaus et al. (2002) */
+ tmp = n[i]*n[k]*CITCOM_DELTA(l,j);
+ tmp += n[j]*n[k]*CITCOM_DELTA(i,l);
+ tmp += n[i]*n[l]*CITCOM_DELTA(k,j);
+ tmp += n[j]*n[l]*CITCOM_DELTA(i,k);
+ tmp /= 2.0;
+ tmp -= 2*n[i]*n[j]*n[k]*n[l];
+ d[i][j][k][l] = tmp;
+ }
+
+}
+
+
+/*
+ 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;
+
+ zero_4x4(c4c);
+
+ for(i1=0;i1<3;i1++)
+ for(i2=0;i2<3;i2++)
+ for(i3=0;i3<3;i3++)
+ for(i4=0;i4<3;i4++)
+ for(j1=0;j1<3;j1++)
+ 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];
+
+}
+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;
+ for(i=0;i<6;i++){
+ for(j=0;j<6;j++)
+ fprintf(out,"%10g ",(fabs(c[i][j])<5e-15)?(0):(c[i][j]));
+ 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];
+}
+
+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.;
+ len += (double)(*x) * (double)(*x);
+ len += (double)(*y) * (double)(*y);
+ len += (double)(*z) * (double)(*z);
+ len = sqrt(len);
+ *x /= len;*y /= len;*z /= len;
+}
+void normalize_vec3d(double *x, double *y, double *z)
+{
+ double len = 0.;
+ len += (*x) * (*x);len += (*y) * (*y);len += (*z) * (*z);
+ 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/CitcomCU/trunk/src/Citcom.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Citcom.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Citcom.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -56,6 +56,7 @@
struct All_variables E;
double time, initial_time, start_time;
+ srand48((long int)-1);
/* parallel_process_initialization(&E,argc,argv); */
E.parallel.me = 0;
@@ -108,7 +109,6 @@
while(E.control.keep_going && (Emergency_stop == 0))
{
process_heating(&E);
-
E.monitor.solution_cycles++;
if(E.monitor.solution_cycles > E.control.print_convergence)
E.control.print_convergence = 1;
@@ -120,20 +120,18 @@
process_temp_field(&E, E.monitor.solution_cycles);
general_stokes_solver(&E);
-
+
if(E.control.composition)
(E.next_buoyancy_field) (&E); /* correct with R-G */
-
/**/ report(&E, "Process results of velocity solver");
process_new_velocity(&E, E.monitor.solution_cycles);
-
if(E.monitor.T_interior > E.monitor.T_interior_max)
{
+ fprintf(stderr, "quit due to maxT = %.4e sub_iteration%d\n", E.monitor.T_interior, E.advection.last_sub_iterations);
fprintf(E.fp, "quit due to maxT = %.4e sub_iteration%d\n", E.monitor.T_interior, E.advection.last_sub_iterations);
parallel_process_termination();
}
-
if(E.parallel.me == 0)
{
fprintf(E.fp, "CPU total = %g & CPU = %g for step %d time = %.4e dt = %.4e maxT = %.4e sub_iteration%d markers=%d\n", CPU_time0() - start_time, CPU_time0() - time, E.monitor.solution_cycles, E.monitor.elapsed_time, E.advection.timestep, E.monitor.T_interior, E.advection.last_sub_iterations, E.advection.markers_g);
@@ -145,8 +143,8 @@
if(E.parallel.me == 0)
{
time = CPU_time0() - initial_time;
- fprintf(E.fp, "Average cpu time taken for velocity step = %f\n", time / ((float)(E.monitor.solution_cycles - 1)));
- fprintf(stderr, "Average cpu time taken for velocity step = %f\n", time / ((float)(E.monitor.solution_cycles - 1)));
+ fprintf(E.fp, "Average cpu time taken for velocity step = %f\n", time / ((float)(E.monitor.solution_cycles)));
+ fprintf(stderr, "Average cpu time taken for velocity step = %f\n", time / ((float)(E.monitor.solution_cycles)));
}
fclose(E.fp);
Modified: mc/3D/CitcomCU/trunk/src/Construct_arrays.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Construct_arrays.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Construct_arrays.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -786,7 +786,9 @@
llayer = layers(E, x3);
if(llayer){
- /* for layers:1-lith,2-upper and 3-lower mantle */
+ /* for layers:1-lith,2-upper and 3-lower mantle
+ or whatever is read in with num_mat and zlayer
+ */
E->mat[el] = llayer;
}
Modified: mc/3D/CitcomCU/trunk/src/Convection.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Convection.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Convection.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -119,6 +119,11 @@
input_int("num_perturbations", &(E->convection.number_of_perturbations), "0,0,32", m);
input_boolean("random_t_init",&E->convection.random_t_init,"off", m);
+#ifndef USE_GGRD
+ if(E->convection.random_t_init)
+ myerror("random T init as per random_t_init only works for GGRD version\n",E);
+#endif
+
input_float_vector("perturbmag", E->convection.number_of_perturbations, E->convection.perturb_mag, m);
input_float_vector("perturbk", E->convection.number_of_perturbations, E->convection.perturb_k, m);
input_float_vector("perturbm", E->convection.number_of_perturbations, E->convection.perturb_mm, m);
@@ -210,9 +215,9 @@
NOTE:
- there's a different routine for GGRD handling in Ggrd_handling
- which has differnet logic for marker init etc. this gets called
- when USE_GGRD is compiled in
+ THERE'S A DIFFERENT ROUTINE FOR GGRD HANDLING IN GGRD_HANDLING
+ WHICH HAS DIFFERNET LOGIC FOR MARKER INIT ETC. THE OTHER ONE GETS
+ CALLED when USE_GGRD is compiled in
=============================== */
Modified: mc/3D/CitcomCU/trunk/src/Drive_solvers.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Drive_solvers.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Drive_solvers.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -40,6 +40,8 @@
#include "element_definitions.h"
#include "global_defs.h"
+int need_to_iterate(struct All_variables *);
+
void general_stokes_solver(struct All_variables *E)
{
//float vmag;
@@ -47,8 +49,7 @@
//double *force, Udot_mag, dUdot_mag;
double Udot_mag, dUdot_mag;
double time;
- //int count, i, j, k;
- int count, i;
+ int i;
static double alpha,alpha1;
static double *oldU;
@@ -59,13 +60,15 @@
//const int dims = E->mesh.nsd;
//const int vpts = vpoints[E->mesh.nsd];
- static int powerlaw;
-
+ int iterate;
+
+ iterate = need_to_iterate(E);
+
if(visits == 0)
{
- powerlaw = (E->viscosity.SDEPV || E->viscosity.BDEPV)?(1):(0);
- if(powerlaw){
+
+ if(iterate){
/* damping factors */
alpha = E->viscosity.sdepv_iter_damp;
alpha1 = 1 - alpha;
@@ -76,7 +79,7 @@
} else{
damp = 0;
}
- } /* end powerlaw */
+ } /* end iterate */
oldU = (double *)malloc(neq * sizeof(double));
for(i = 0; i < neq; i++)
oldU[i] = 0.0;
@@ -96,7 +99,6 @@
time = CPU_time0();
velocities_conform_bcs(E, E->U);
-
assemble_forces(E, 0);
/*
@@ -107,19 +109,18 @@
}
*/
- count = 1;
+ E->monitor.visc_iter_count = 0;
do
{
-
if(E->viscosity.update_allowed)
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(powerlaw)
+ E->monitor.visc_iter_count++;
+
+ if(iterate)
{
if(damp){
@@ -140,15 +141,32 @@
if(E->control.sdepv_print_convergence && (E->parallel.me == 0))
{
- fprintf(stderr, "Stress dependent viscosity: DUdot = %.4e (%.4e) for iteration %d\n", dUdot_mag, Udot_mag, count);
- fprintf(E->fp, "Stress dependent viscosity: DUdot = %.4e (%.4e) for iteration %d\n", dUdot_mag, Udot_mag, count);
+ fprintf(stderr, "Stress dependent viscosity: DUdot = %.4e (%.4e) for iteration %d\n", dUdot_mag, Udot_mag, E->monitor.visc_iter_count);
+ fprintf(E->fp, "Stress dependent viscosity: DUdot = %.4e (%.4e) for iteration %d\n", dUdot_mag, Udot_mag, E->monitor.visc_iter_count);
fflush(E->fp);
}
- count++;
+ E->monitor.visc_iter_count++;
} /* end for SDEPV / BDEPV */
- } while((count < 50) && (dUdot_mag > E->viscosity.sdepv_misfit) && powerlaw);
+ } while((E->monitor.visc_iter_count < 50) && (dUdot_mag > E->viscosity.sdepv_misfit) && iterate);
free((void *)delta_U);
return;
}
+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.BDEPV)?(TRUE):(FALSE);
+ }
+#else
+ /* regular operation */
+ return ((E->viscosity.SDEPV || E->viscosity.BDEPV)?(TRUE):(FALSE));
+#endif
+}
Modified: mc/3D/CitcomCU/trunk/src/Element_calculations.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Element_calculations.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Element_calculations.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -45,6 +45,11 @@
#include <sys/resource.h>
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+#include "anisotropic_viscosity.h"
+#endif
+
+
/* *INDENT-OFF* */
static double Dl[5][5] = { {0.0, 0.0, 0.0, 0.0, 0.0},
{0.0, 1.0, 1.0, 0.0, 1.0},
@@ -80,7 +85,6 @@
//const int npno = E->lmesh.npno;
const int nel = E->lmesh.nel;
const int lev = E->mesh.levmax;
-
for(a = 0; a < neq; a++)
E->F[a] = 0.0;
@@ -127,160 +131,248 @@
void get_elt_k(struct All_variables *E, int el, double elt_k[24 * 24], int lev, int iconv)
{
- double bdbmu[4][4];
- //double bdbl[4][4];
+ double bdbmu[4][4];
- double rtf[4][9], W[9], ra[9], si[9], ct[9];
- //struct Shape_function GN;
- //struct Shape_function_dA dOmega;
- //struct Shape_function_dx GNx;
- static struct CC Cc;
- static struct CCX Ccx;
+ double rtf[4][9], W[9];
+ static struct CC Cc;
+ static struct CCX Ccx;
- //int p1[9], pn, qn, ad, bd;
- int pn, qn, ad, bd;
- //int nodea, nodeb, a, b, i, j, k, p, q, nint;
- int a, b, i, j, k;
- //double RM2[9], RMP[9], r;
- //double Visc, visc[9], temp;
- double temp;
- float gnx0, gnx1, gnx2;
- double shp, cc1, cc2, cc3;
+ int pn, qn, ad, bd;
+ int a, b, i, j, k;
+ double temp;
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ double D[VPOINTS3D+1][6][6],btmp[6];
+ int l1,l2;
+#endif
+ int off;
+ double ba[9][4][9][7]; /* flipped around for improved speed (different from CitcomS) */
- double ba[7][4][9][9];
+ const int n = loc_mat_size[E->mesh.nsd];
+ const int vpts = vpoints[E->mesh.nsd];
+ //const int ppts = ppoints[E->mesh.nsd];
+ const int ends = enodes[E->mesh.nsd];
+ const int dims = E->mesh.nsd;
+ //const int dofs = E->mesh.dof;
+ //const int sphere_key = 1;
+ //const double zero = 0.0;
+ const double one = 1.0;
+ const double two = 2.0;
- const int n = loc_mat_size[E->mesh.nsd];
- const int vpts = vpoints[E->mesh.nsd];
- //const int ppts = ppoints[E->mesh.nsd];
- const int ends = enodes[E->mesh.nsd];
- const int dims = E->mesh.nsd;
- //const int dofs = E->mesh.dof;
- //const int sphere_key = 1;
- //const double zero = 0.0;
- const double one = 1.0;
- const double two = 2.0;
- for(k = 1; k <= vpts; k++)
- {
- W[k] = g_point[k].weight[dims - 1] * E->GDA[lev][el].vpt[k] * E->EVI[lev][(el - 1) * vpts + k];
- }
+ if(E->control.Rsphere) /* need rtf for spherical */
+ get_rtf(E, el, 0, rtf, lev);
+ for(k = 1; k <= vpts; k++){
+ off = (el-1)*vpts+k;
+ W[k] = g_point[k].weight[dims - 1] * E->GDA[lev][el].vpt[k] * E->EVI[lev][off];
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity){
+ /* allow for a possibly anisotropic viscosity, only used if switched on */
+ get_constitutive(D[k],lev,off,rtf[1][k],rtf[2][k],(E->control.Rsphere),E);
+ }
+#endif
+ }
- if(E->control.Rsphere)
- {
- get_rtf(E, el, 0, rtf, lev);
+ if(E->control.Rsphere){
+ if(iconv == 1 || ((iconv == 0) && (el - 1) % E->lmesh.ELZ[lev] == 0))
+ construct_c3x3matrix_el(E, el, &Cc, &Ccx, lev, 0);
+ get_ba(&(E->N),&(E->GNX[lev][el]),&Cc, &Ccx,rtf, E->mesh.nsd,TRUE,ba);
+ }else{ /* end for sphere */
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity){ /* only anisotropic cartesian uses ba[] */
+ if(iconv == 1 || ((iconv == 0) && (el - 1) % E->lmesh.ELZ[lev] == 0))
+ get_ba(&(E->N),&(E->GNX[lev][el]),&Cc, &Ccx,rtf, E->mesh.nsd, FALSE,ba);
+ }
+#endif
+ ; /* only anisotropic cartesian needs ba factors */
+ }
- if(iconv == 1 || ((iconv == 0) && (el - 1) % E->lmesh.ELZ[lev] == 0))
- construct_c3x3matrix_el(E, el, &Cc, &Ccx, lev, 0);
+ for(a = 1; a <= ends; a++)
+ for(b = a; b <= ends; b++){
+ bdbmu[1][1] = bdbmu[1][2] = bdbmu[1][3] =
+ bdbmu[2][1] = bdbmu[2][2] = bdbmu[2][3] =
+ bdbmu[3][1] = bdbmu[3][2] = bdbmu[3][3] = 0.0;
+#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 <= VPOINTS3D; k++){
+ /* note that D is in 0,...,N-1 convention */
+ for(l1=0;l1 < 6;l1++){ /* compute D*B */
+ btmp[l1] = ( D[k][l1][0] * ba[b][j][k][1] + D[k][l1][1] * ba[b][j][k][2] + D[k][l1][2] * ba[b][j][k][3] +
+ D[k][l1][3] * ba[b][j][k][4] + D[k][l1][4] * ba[b][j][k][5] + D[k][l1][5] * ba[b][j][k][6] );
+ }
+ /* compute B^T (D*B) */
+ bdbmu[i][j] += W[k] * ( ba[a][i][k][1]*btmp[0] + ba[a][i][k][2]*btmp[1] + ba[a][i][k][3]*btmp[2]+
+ ba[a][i][k][4]*btmp[3] + ba[a][i][k][5]*btmp[4] + ba[a][i][k][6]*btmp[5]);
+ }
+ }else{
+#endif
+ /*
- for(k = 1; k <= vpts; k++)
- {
- ra[k] = rtf[3][k];
- si[k] = one / sin(rtf[1][k]);
- ct[k] = cos(rtf[1][k]) * si[k];
- }
+ isotropic branch, unchanged from previous version where only
+ the spherical version uses B
- for(a = 1; a <= ends; a++)
- for(i = 1; i <= dims; i++)
- for(k = 1; k <= VPOINTS3D; k++)
- {
- gnx0 = E->GNX[lev][el].vpt[GNVXINDEX(0, a, k)];
- gnx1 = E->GNX[lev][el].vpt[GNVXINDEX(1, a, k)];
- gnx2 = E->GNX[lev][el].vpt[GNVXINDEX(2, a, k)];
- shp = E->N.vpt[GNVINDEX(a, k)];
- cc1 = Cc.vpt[BVINDEX(1, i, a, k)];
- cc2 = Cc.vpt[BVINDEX(2, i, a, k)];
- cc3 = Cc.vpt[BVINDEX(3, i, a, k)];
+ */
+ if(E->control.CART3D){
+ /* cartesian isotropic does not use ba[] */
+ for(k = 1; k <= VPOINTS3D; k++){
+ bdbmu[1][1] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(0, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(0, b, k)];
+ bdbmu[1][2] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(1, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(0, b, k)];
+ bdbmu[1][3] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(2, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(0, b, k)];
+ bdbmu[2][1] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(0, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(1, b, k)];
+ bdbmu[2][2] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(1, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(1, b, k)];
+ bdbmu[2][3] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(2, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(1, b, k)];
+ bdbmu[3][1] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(0, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(2, b, k)];
+ bdbmu[3][2] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(1, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(2, b, k)];
+ bdbmu[3][3] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(2, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(2, b, k)];
+ }
+
+ temp = 0.0;
+ for(k = 1; k <= VPOINTS3D; k++){
+ temp += W[k] * (E->GNX[lev][el].vpt[GNVXINDEX(0, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(0, b, k)] +
+ E->GNX[lev][el].vpt[GNVXINDEX(1, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(1, b, k)] +
+ E->GNX[lev][el].vpt[GNVXINDEX(2, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(2, b, k)]);
+ }
+
+ bdbmu[1][1] += temp;
+ bdbmu[2][2] += temp;
+ bdbmu[3][3] += temp;
+
+ /* end for Cart3D */
+ }else if(E->control.Rsphere){
+ for(i = 1; i <= dims; i++)
+ for(j = 1; j <= dims; j++)
+ for(k = 1; k <= VPOINTS3D; k++){
+ bdbmu[i][j] += W[k] * (two * (ba[a][i][k][1] * ba[b][j][k][1] +
+ ba[a][i][k][2] * ba[b][j][k][2] +
+ ba[a][i][k][3] * ba[b][j][k][3]) +
+ ba[a][i][k][4] * ba[b][j][k][4] +
+ ba[a][i][k][5] * ba[b][j][k][5] +
+ ba[a][i][k][6] * ba[b][j][k][6]);
+ }
+ } /* end for Sphere */
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ } /* end isotropic */
+#endif
+
+ ad = dims * (a - 1);
+ bd = dims * (b - 1);
+ pn = ad * n + bd;
+ qn = bd * n + ad;
+
+ elt_k[pn] = bdbmu[1][1]; /* above */
+ elt_k[pn + 1] = bdbmu[1][2];
+ elt_k[pn + 2] = bdbmu[1][3];
+ elt_k[pn + n] = bdbmu[2][1];
+ elt_k[pn + n + 1] = bdbmu[2][2];
+ elt_k[pn + n + 2] = bdbmu[2][3];
+ elt_k[pn + 2 * n] = bdbmu[3][1];
+ elt_k[pn + 2 * n + 1] = bdbmu[3][2];
+ elt_k[pn + 2 * n + 2] = bdbmu[3][3];
+
+ elt_k[qn] = bdbmu[1][1]; /* below diag */
+ elt_k[qn + 1] = bdbmu[2][1];
+ elt_k[qn + 2] = bdbmu[3][1];
+ elt_k[qn + n] = bdbmu[1][2];
+ elt_k[qn + n + 1] = bdbmu[2][2];
+ elt_k[qn + n + 2] = bdbmu[3][2];
+ elt_k[qn + 2 * n] = bdbmu[1][3];
+ elt_k[qn + 2 * n + 1] = bdbmu[2][3];
+ elt_k[qn + 2 * n + 2] = bdbmu[3][3];
+ /**/
+ } /* Sum over all the a,b's to obtain full elt_k matrix */
+
+ return;
+}
+/*
- ba[1][i][a][k] = (gnx0 * cc1 + shp * Ccx.vpt[BVXINDEX(1, i, 1, a, k)] + shp * cc3) * ra[k];
+compute the displacement gradient matrix B
- ba[2][i][a][k] = (shp * cc1 * ct[k] + shp * cc3 + (gnx1 * cc2 + shp * Ccx.vpt[BVXINDEX(2, i, 2, a, k)]) * si[k]) * ra[k];
+*/
+void get_ba(struct Shape_function *N,struct Shape_function_dx *GNx,
+ struct CC *cc, struct CCX *ccx, double rtf[4][9],
+ int dims,int spherical, double ba[9][4][9][7])
+{
- ba[3][i][a][k] = gnx2 * cc3;
+ int i,k,a,vpts,ends;
+ double ra[9], si[9], ct[9];
+ const double one = 1.0;
+ const double two = 2.0;
+ float gnx0, gnx1, gnx2;
+ double shp, cc1, cc2, cc3;
+
+ vpts = vpoints[dims];
+ ends = enodes[dims];
- ba[4][i][a][k] = (gnx0 * cc2 + shp * Ccx.vpt[BVXINDEX(2, i, 1, a, k)] - shp * cc2 * ct[k] + (gnx1 * cc1 + shp * Ccx.vpt[BVXINDEX(1, i, 2, a, k)]) * si[k]) * ra[k];
- ba[5][i][a][k] = gnx2 * cc1 + (gnx0 * cc3 + shp * (Ccx.vpt[BVXINDEX(3, i, 1, a, k)] - cc1)) * ra[k];
+ if(spherical){
+ for(k = 1; k <= vpts; k++){
+ ra[k] = rtf[3][k];
+ si[k] = one / sin(rtf[1][k]);
+ ct[k] = cos(rtf[1][k]) * si[k];
+ }
+ for(a = 1; a <= ends; a++)
+ for(i = 1; i <= dims; i++)
+ for(k = 1; k <= VPOINTS3D; k++){
+ gnx0 = GNx->vpt[GNVXINDEX(0, a, k)];
+ gnx1 = GNx->vpt[GNVXINDEX(1, a, k)];
+ gnx2 = GNx->vpt[GNVXINDEX(2, a, k)];
+ shp = N->vpt[GNVINDEX(a, k)];
+ cc1 = cc->vpt[BVINDEX(1, i, a, k)];
+ cc2 = cc->vpt[BVINDEX(2, i, a, k)];
+ cc3 = cc->vpt[BVINDEX(3, i, a, k)];
+
+ ba[a][i][k][1] = (gnx0 * cc1 + shp * ccx->vpt[BVXINDEX(1, i, 1, a, k)] + shp * cc3) * ra[k];
+
+ ba[a][i][k][2] = (shp * cc1 * ct[k] + shp * cc3 + (gnx1 * cc2 + shp * ccx->vpt[BVXINDEX(2, i, 2, a, k)]) * si[k]) * ra[k];
+
+ ba[a][i][k][3] = gnx2 * cc3;
+
+ ba[a][i][k][4] = (gnx0 * cc2 + shp * ccx->vpt[BVXINDEX(2, i, 1, a, k)] - shp * cc2 * ct[k] + (gnx1 * cc1 + shp * ccx->vpt[BVXINDEX(1, i, 2, a, k)]) * si[k]) * ra[k];
+
+ ba[a][i][k][5] = gnx2 * cc1 + (gnx0 * cc3 + shp * (ccx->vpt[BVXINDEX(3, i, 1, a, k)] - cc1)) * ra[k];
+
+ ba[a][i][k][6] = gnx2 * cc2 - ra[k] * shp * cc2 + (gnx1 * cc3 + shp * ccx->vpt[BVXINDEX(3, i, 2, a, k)]) * si[k] * ra[k];
+ }
+ }else{ /* cartesian */
+ for(a = 1; a <= ends; a++)
+ for(k = 1; k <= VPOINTS3D; k++){
+ gnx0 = GNx->vpt[GNVXINDEX(0, a, k)];
+ gnx1 = GNx->vpt[GNVXINDEX(1, a, k)];
+ gnx2 = GNx->vpt[GNVXINDEX(2, a, k)];
- ba[6][i][a][k] = gnx2 * cc2 - ra[k] * shp * cc2 + (gnx1 * cc3 + shp * Ccx.vpt[BVXINDEX(3, i, 2, a, k)]) * si[k] * ra[k];
- }
+ ba[a][1][k][1] = gnx0;
+ ba[a][2][k][1] = 0.;
+ ba[a][3][k][1] = 0.;
+ ba[a][1][k][2] = 0.;
+ ba[a][2][k][2] = gnx1;
+ ba[a][3][k][2] = 0.;
- } /* end for sphere */
+ ba[a][1][k][3] = 0.;
+ ba[a][2][k][3] = 0.;
+ ba[a][3][k][3] = gnx2;
+ ba[a][1][k][4] = gnx1;
+ ba[a][2][k][4] = gnx0;
+ ba[a][3][k][4] = 0.;
- for(a = 1; a <= ends; a++)
- for(b = a; b <= ends; b++)
- {
- bdbmu[1][1] = bdbmu[1][2] = bdbmu[1][3] = bdbmu[2][1] = bdbmu[2][2] = bdbmu[2][3] = bdbmu[3][1] = bdbmu[3][2] = bdbmu[3][3] = 0.0;
+ ba[a][1][k][5] = gnx2;
+ ba[a][2][k][5] = 0.;
+ ba[a][3][k][5] = gnx0;
- if(E->control.CART3D)
- {
- for(k = 1; k <= VPOINTS3D; k++)
- {
- bdbmu[1][1] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(0, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(0, b, k)];
- bdbmu[1][2] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(1, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(0, b, k)];
- bdbmu[1][3] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(2, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(0, b, k)];
- bdbmu[2][1] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(0, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(1, b, k)];
- bdbmu[2][2] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(1, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(1, b, k)];
- bdbmu[2][3] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(2, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(1, b, k)];
- bdbmu[3][1] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(0, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(2, b, k)];
- bdbmu[3][2] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(1, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(2, b, k)];
- bdbmu[3][3] += W[k] * E->GNX[lev][el].vpt[GNVXINDEX(2, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(2, b, k)];
- }
+ ba[a][1][k][6] = 0.;
+ ba[a][2][k][6] = gnx2;
+ ba[a][3][k][6] = gnx1;
- temp = 0.0;
- for(k = 1; k <= VPOINTS3D; k++)
- temp += W[k] * (E->GNX[lev][el].vpt[GNVXINDEX(0, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(0, b, k)] + E->GNX[lev][el].vpt[GNVXINDEX(1, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(1, b, k)] + E->GNX[lev][el].vpt[GNVXINDEX(2, a, k)] * E->GNX[lev][el].vpt[GNVXINDEX(2, b, k)]);
+ }
+ } /* end Caretsian */
- bdbmu[1][1] += temp;
- bdbmu[2][2] += temp;
- bdbmu[3][3] += temp;
-
- } /* end for Cart3D */
-
- else if(E->control.Rsphere)
- {
- for(i = 1; i <= dims; i++)
- for(j = 1; j <= dims; j++)
- for(k = 1; k <= VPOINTS3D; k++)
- {
- bdbmu[i][j] += W[k] * (two * (ba[1][i][a][k] * ba[1][j][b][k] + ba[2][i][a][k] * ba[2][j][b][k] + ba[3][i][a][k] * ba[3][j][b][k]) + ba[4][i][a][k] * ba[4][j][b][k] + ba[5][i][a][k] * ba[5][j][b][k] + ba[6][i][a][k] * ba[6][j][b][k]);
- }
- } /* end for Sphere */
-
- ad = dims * (a - 1);
- bd = dims * (b - 1);
- pn = ad * n + bd;
- qn = bd * n + ad;
-
- elt_k[pn] = bdbmu[1][1]; /* above */
- elt_k[pn + 1] = bdbmu[1][2];
- elt_k[pn + 2] = bdbmu[1][3];
- elt_k[pn + n] = bdbmu[2][1];
- elt_k[pn + n + 1] = bdbmu[2][2];
- elt_k[pn + n + 2] = bdbmu[2][3];
- elt_k[pn + 2 * n] = bdbmu[3][1];
- elt_k[pn + 2 * n + 1] = bdbmu[3][2];
- elt_k[pn + 2 * n + 2] = bdbmu[3][3];
-
- elt_k[qn] = bdbmu[1][1]; /* below diag */
- elt_k[qn + 1] = bdbmu[2][1];
- elt_k[qn + 2] = bdbmu[3][1];
- elt_k[qn + n] = bdbmu[1][2];
- elt_k[qn + n + 1] = bdbmu[2][2];
- elt_k[qn + n + 2] = bdbmu[3][2];
- elt_k[qn + 2 * n] = bdbmu[1][3];
- elt_k[qn + 2 * n + 1] = bdbmu[2][3];
- elt_k[qn + 2 * n + 2] = bdbmu[3][3];
- /**/} /* Sum over all the a,b's to obtain full elt_k matrix */
-
- return;
}
-
/* =============================================
* General calling function for del_squared:
* according to whether it should be element by
Modified: mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -56,7 +56,7 @@
unsigned int type;
FILE *fp;
- double x1, y1, z1, con;
+ double x1, y1, z1, con, top_t, bot_t;
int noz2, nfz, in1, in2, in3, instance, nox, noy, noz,nxs,setflag;
char input_s[200], output_file[255];
@@ -67,7 +67,7 @@
/* twb additions */
double rho_prem;
char pfile[1000];
- double t1,f1,r1,tgrad,tbot,tadd,tz,tmean;
+ double t1,f1,r1,tgrad,tadd,tz,tmean;
/* for dealing with several processors */
MPI_Status mpi_stat;
@@ -84,6 +84,10 @@
setflag = 0;
+ /* top and bottom temperatures for initial assign (only use if
+ temperatures are set, else defaults */
+ bot_t = (E->mesh.bottbc) ? E->control.TBCbotval : 1.0;
+ top_t = (E->mesh.toptbc) ? E->control.TBCtopval : 1.0;
for(i=1;i<=E->lmesh.nno;i++) {
/* init as zeros */
@@ -161,26 +165,22 @@
interpolate densities to temperature given PREM variations
*/
- if(E->mesh.bottbc != 0){
- /* bottom has specified temperature */
- tbot = E->control.TBCbotval;
- }else{
- /*
- bottom has specified heat flux
- start with unity bottom temperature
- */
- tbot = 1.0;
- }
- tmean = (tbot + E->control.TBCtopval)/2.0 + E->control.ggrd.temp.offset;
+
+ tmean = (top_t + bot_t)/2.0 + E->control.ggrd.temp.offset;
if(E->parallel.me == 0)
fprintf(stderr,"tinit: top: %g bot: %g mean: %g scale: %g PREM %i\n",
- E->control.TBCtopval,tbot,tmean,E->control.ggrd.temp.scale,
+ top_t,bot_t,tmean,E->control.ggrd.temp.scale,
E->control.ggrd.temp.scale_with_prem);
for(i=1;i<=noy;i++)
for(j=1;j<=nox;j++)
for(k=1;k<=noz;k++) {
node=k+(j-1)*noz+(i-1)*nox*noz; /* offset */
if(E->control.slab_slice){
+ /*
+
+ slab slice input
+
+ */
if(E->control.Rsphere) {
if(E->SX[1][node] <= E->control.slab_theta_bound)
/* spherical interpolation */
@@ -210,7 +210,7 @@
}else{
/*
- 3-D
+ 3-D grid model input
*/
if(E->control.Rsphere) /* spherical interpolation */
@@ -284,7 +284,9 @@
add a linear profile and potentially perturbations to the temperature
+
*/
+
if(E->control.CART3D)
{
@@ -296,15 +298,9 @@
x1 = E->X[1][node];
z1 = E->X[3][node];
y1 = E->X[2][node];
- /* linear gradient */
- if(E->mesh.bottbc){
- tz = z1* (E->control.TBCtopval - E->control.TBCbotval) +
- E->control.TBCbotval;
- }else{
- tz = z1* (E->control.TBCtopval - 1) + 1;
- }
+ /* linear gradient assuming unity box height */
+ tz = z1 * (top_t - bot_t) + bot_t;
E->T[node] += tz;
-
if(E->convection.random_t_init){
/* random init */
E->T[node] += (-0.5 + drand48()) * E->convection.perturb_mag[0];
@@ -338,8 +334,8 @@
z1 = E->SX[3][node];
y1 = E->SX[2][node];
/* linear */
- tz = (z1 - E->sphere.ri)/(E->sphere.ro - E->sphere.ri) *
- (E->control.TBCtopval-E->control.TBCbotval) + E->control.TBCbotval;
+ tz = (z1 - E->sphere.ri)/(E->sphere.ro - E->sphere.ri) * (top_t-bot_t) + bot_t;
+
E->T[node] += tz;
if(E->convection.random_t_init){
/* random init */
@@ -492,3 +488,253 @@
return;
}
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+/*
+
+
+read in anisotropic viscosity from a directory which holds
+
+
+vis2.grd for the viscosity factors, read in log10(eta_S/eta)
+
+nr.grd, nt.grd, np.grd for the directors
+
+*/
+void ggrd_read_anivisc_from_file(struct All_variables *E)
+{
+ 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 llayer,nox,noy,noz,level,lselect,idim,elx,ely,elz;
+ char gmt_string[10],char_dummy;
+ double vis2,log_vis,ntheta,nphi,nr;
+ float rout[3],xloc[3];
+ double cvec[3];
+ float base[9];
+ char tfilename[1000];
+ static ggrd_boolean shift_to_pos_lon = FALSE;
+ const int dims=E->mesh.nsd;
+ 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];
+
+
+ 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);
+
+ /* isotropic default */
+ for(i=E->mesh.levmin;i <= E->mesh.levmax;i++){
+ nel = E->lmesh.NEL[i];
+ for(k=1;k <= nel;k++){
+ for(l=1;l <= vpts;l++){ /* assign to all integration points */
+ ind = (k-1)*vpts + l;
+ E->EVI2[i][ind] = 0.0;
+ E->EVIn1[i][ind] = 1.0; E->EVIn2[i][ind] = E->EVIn3[i][ind] = 0.0;
+ }
+ }
+ }
+ sprintf(gmt_string,""); /* regional */
+
+ /*
+
+ initialization steps
+
+ */
+ if(E->parallel.me > 0) /* wait for previous processor */
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1),
+ 0, MPI_COMM_WORLD, &mpi_stat);
+ /*
+ read in the material file(s)
+ */
+ vis2_grd = (struct ggrd_gt *)calloc(1,sizeof(struct ggrd_gt));
+ nphi_grd = (struct ggrd_gt *)calloc(1,sizeof(struct ggrd_gt));
+ nr_grd = (struct ggrd_gt *)calloc(1,sizeof(struct ggrd_gt));
+ ntheta_grd = (struct ggrd_gt *)calloc(1,sizeof(struct ggrd_gt));
+
+
+ if(E->parallel.me==0)
+ 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->monitor.length_scale*E->viscosity.zbase_layer[E->viscosity.anivisc_layer - 1],
+ ("regional"));
+ else
+ fprintf(stderr,"ggrd_read_anivisc_from_file: initializing, assigning to all elements between %g and %g km, input is %s\n",
+ E->monitor.length_scale*((E->viscosity.anivisc_layer<-1)?(E->viscosity.zbase_layer[-E->viscosity.anivisc_layer - 2]):(0)),
+ E->monitor.length_scale*E->viscosity.zbase_layer[-E->viscosity.anivisc_layer - 1],
+ ("regional"));
+
+ /*
+ read viscosity ratio, and east/north direction of normal azimuth
+ */
+ /* viscosity factor */
+ sprintf(tfilename,"%s/vis2.grd",E->viscosity.anisotropic_init_dir);
+ if(ggrd_grdtrack_init_general(FALSE,tfilename,"",gmt_string,
+ vis2_grd,(E->parallel.me == 0),FALSE))
+ myerror("ggrd init error",E);
+ /* n_r */
+ if(E->control.CART3D)
+ sprintf(tfilename,"%s/nz.grd",E->viscosity.anisotropic_init_dir);
+ else
+ sprintf(tfilename,"%s/nr.grd",E->viscosity.anisotropic_init_dir);
+
+ if(ggrd_grdtrack_init_general(FALSE,tfilename,"",gmt_string,
+ nr_grd,(E->parallel.me == 0),FALSE))
+ myerror("ggrd init error",E);
+ if(E->control.CART3D)
+ /* n_theta */
+ sprintf(tfilename,"%s/ny.grd",E->viscosity.anisotropic_init_dir);
+ else
+ sprintf(tfilename,"%s/nt.grd",E->viscosity.anisotropic_init_dir);
+
+ if(ggrd_grdtrack_init_general(FALSE,tfilename,"",gmt_string,
+ ntheta_grd,(E->parallel.me == 0),FALSE))
+ myerror("ggrd init error",E);
+ if(E->control.CART3D)
+ /* n_phi */
+ sprintf(tfilename,"%s/nx.grd",E->viscosity.anisotropic_init_dir);
+ else
+ sprintf(tfilename,"%s/np.grd",E->viscosity.anisotropic_init_dir);
+
+ if(ggrd_grdtrack_init_general(FALSE,tfilename,"",gmt_string,
+ nphi_grd,(E->parallel.me == 0),FALSE))
+ myerror("ggrd init error",E);
+
+ /* done */
+ if(E->parallel.me < E->parallel.nproc-1){ /* tell the next proc to go ahead */
+ mpi_rc = MPI_Send(&mpi_success_message, 1,
+ MPI_INT, (E->parallel.me+1), 0, MPI_COMM_WORLD);
+ }else{
+ fprintf(stderr,"ggrd_read_anivisc_from_file: last processor done with ggrd init\n");
+ fprintf(stderr,"ggrd_read_anivisc_from_file: WARNING: assuming a regular grid geometry\n");
+ }
+ /*
+
+ loop through all elements and assign
+
+ */
+ for (j=1;j <= elz;j++) { /* this assumes a regular grid sorted as in (1)!!! */
+ if(((E->viscosity.anivisc_layer > 0)&&(E->mat[j] <= E->viscosity.anivisc_layer))||
+ ((E->viscosity.anivisc_layer < 0)&&(E->mat[j] == -E->viscosity.anivisc_layer))){
+ /* within top layers */
+ for (k=1;k <= ely;k++){
+ for (i=1;i <= elx;i++) {
+ /* eq.(1) */
+ el = j + (i-1) * elz + (k-1)*elxlz;
+ /*
+ find average coordinates
+ */
+ if(E->control.CART3D){
+ rout[0]=rout[1]=rout[2]=0.0;
+ for(inode=1;inode <= 4;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.;
+ 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]);
+ parallel_process_termination();
+ }
+ /* read in cartesian vector */
+ if(!ggrd_grdtrack_interpolate_xy((double)rout[0],(double)rout[1],nphi_grd,cvec,FALSE)){
+ fprintf(stderr,"ggrd_read_anivisc_from_file: interpolation error at x: %g y: %g\n",
+ rout[0],rout[1]);
+ parallel_process_termination();
+ }
+ if(!ggrd_grdtrack_interpolate_xy((double)rout[0],(double)rout[1],ntheta_grd,(cvec+1),FALSE)){
+ fprintf(stderr,"ggrd_read_anivisc_from_file: interpolation error at x: %g y: %g\n",
+ rout[0],rout[1]);
+ parallel_process_termination();
+ }
+ if(!ggrd_grdtrack_interpolate_xy((double)rout[0],(double)rout[1],nr_grd,(cvec+2),FALSE)){
+ fprintf(stderr,"ggrd_read_anivisc_from_file: interpolation error at x: %g y: %g\n",
+ rout[0],rout[1]);
+ parallel_process_termination();
+ }
+ normalize_vec3d(cvec,(cvec+1),(cvec+2));
+ }else{
+ /* spherical */
+ xloc[0] = xloc[1] = xloc[2] = 0.0;
+ for(inode=1;inode <= 4;inode++){
+ ind = E->ien[el].node[inode];
+ rtp2xyz((float)E->SX[3][ind],(float)E->SX[1][ind],
+ (float)E->SX[2][ind],rout);
+ xloc[0] += rout[0];
+ xloc[1] += rout[1];
+ xloc[2] += rout[2];
+ }
+ xloc[0]/=4.;xloc[1]/=4.;xloc[2]/=4.;
+ xyz2rtp(xloc[0],xloc[1],xloc[2],rout); /* convert to spherical */
+
+ /* vis2 */
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],
+ 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();
+ }
+ /* nr */
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],
+ nr_grd,&nr,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();
+ }
+ /* ntheta */
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],
+ ntheta_grd,&ntheta,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();
+ }
+ /* nphi */
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],
+ nphi_grd,&nphi,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();
+ }
+ normalize_vec3d(&nr,&ntheta,&nphi);
+ calc_cbase_at_tp(rout[1],rout[2],base); /* convert from
+ spherical
+ coordinates to
+ Cartesian */
+ convert_pvec_to_cvec(nr,ntheta,nphi,base,xloc);
+ cvec[0]=xloc[0];cvec[1]=xloc[1];cvec[2]=xloc[2];
+ }
+ /* 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.levmax][ind] = vis2;
+ E->EVIn1[E->mesh.levmax][ind] = cvec[0];
+ E->EVIn2[E->mesh.levmax][ind] = cvec[1];
+ E->EVIn3[E->mesh.levmax][ind] = cvec[2];
+ }
+ }
+ }
+ } /* end insize lith */
+ } /* end elz loop */
+
+
+ ggrd_grdtrack_free_gstruc(vis2_grd);
+ ggrd_grdtrack_free_gstruc(nr_grd);
+ ggrd_grdtrack_free_gstruc(ntheta_grd);
+ ggrd_grdtrack_free_gstruc(nphi_grd);
+
+
+}
+
+
+#endif /* for ANISOTROPIC */
+
Modified: mc/3D/CitcomCU/trunk/src/Instructions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Instructions.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Instructions.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -172,7 +172,7 @@
void allocate_common_vars(struct All_variables *E)
{
//int nox, noy, noz, i, j, l, nno_l, npno_l, nozl, nnov_l, nxyz;
- int nox, noy, noz, i, j, l, nxyz;
+ int nox, noy, noz, i, j, l, nxyz,nel,nno;
E->mesh.fnodal_malloc_size = (E->lmesh.nno + 2) * sizeof(float);
E->mesh.dnodal_malloc_size = (E->lmesh.nno + 2) * sizeof(double);
@@ -282,7 +282,38 @@
E->BPI[i] = (double *)malloc((E->lmesh.NPNO[i] + 1) * sizeof(double));
E->control.B_is_good[i] = 0;
}
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity){ /* any anisotropic
+ viscosity */
+ for(i=E->mesh.levmin;i<=E->mesh.levmax;i++){
+ nel = E->lmesh.NEL[i];
+ nno = E->lmesh.NNO[i];
+ E->EVI2[i] = (float *) malloc((nel+1)*vpoints[E->mesh.nsd]*sizeof(float));
+ E->EVIn1[i] = (float *) malloc((nel+1)*vpoints[E->mesh.nsd]*sizeof(float));
+ E->EVIn2[i] = (float *) malloc((nel+1)*vpoints[E->mesh.nsd]*sizeof(float));
+ E->EVIn3[i] = (float *) malloc((nel+1)*vpoints[E->mesh.nsd]*sizeof(float));
+
+ E->VI2[i] = (float *) malloc((nno+1)*sizeof(float));
+ E->VIn1[i] = (float *) malloc((nno+1)*sizeof(float));
+ E->VIn2[i] = (float *) malloc((nno+1)*sizeof(float));
+ E->VIn3[i] = (float *) malloc((nno+1)*sizeof(float));
+ if((!(E->EVI2[i]))||(!(E->VI2[i]))||
+ (!(E->EVIn1[i]))||(!(E->EVIn2[i]))||(!(E->EVIn3[i]))||
+ (!(E->VIn1[i]))||(!(E->VIn2[i]))||(!(E->VIn3[i]))){
+ fprintf(stderr, "Error: Cannot allocate anisotropic visc memory, rank=%d\n",
+ E->parallel.me);
+ parallel_process_termination();
+ }
+ }
+ E->viscosity.anisotropic_viscosity_init = FALSE;
+ if(E->parallel.me == 0)
+ fprintf(stderr,"allocated for anisotropic viscosity (%s) levmax %i\n",
+ (E->viscosity.allow_anisotropic_viscosity == 1)?("orthotropic"):("transversely isotropic"),
+ E->mesh.levmax);
+ }
+#endif
+
E->temp = (double *)malloc((E->lmesh.NEQ[E->mesh.levmax] + 2) * sizeof(double));
E->Element = (unsigned int *)malloc((E->lmesh.nel + 2) * sizeof(unsigned int));
@@ -703,6 +734,8 @@
set_2dc_defaults(E);
}
+
+
input_string("Solver", E->control.SOLVER_TYPE, NULL, m);
if(strcmp(E->control.SOLVER_TYPE, "cgrad") == 0)
{
@@ -789,25 +822,8 @@
input_int("obs_maxlongk", &(E->slice.maxlong), "100,1", m);
input_int("obs_minlongk", &(E->slice.minlong), "1,1", m);
- /* for layers */
- E->viscosity.zlm = 1.0;
- E->viscosity.z410 = 1.0;
- E->viscosity.zlith = 0.0;
- if(E->control.CART3D)
- {
- input_float("z_lmantle", &(E->viscosity.zlm), "1.0", m);
- input_float("z_410", &(E->viscosity.z410), "1.0", m);
- input_float("z_lith", &(E->viscosity.zlith), "0.0", m);
- }
- else if(E->control.Rsphere)
- {
- input_float("r_lmantle", &(E->viscosity.zlm), "1.0", m);
- input_float("r_410", &(E->viscosity.z410), "1.0", m);
- input_float("r_lith", &(E->viscosity.zlith), "0.0", m);
- }
-
#ifdef USE_GGRD
/* ggrd control */
ggrd_init_master(&(E->control.ggrd));
@@ -986,8 +1002,10 @@
E->monitor.tau_scale = (E->data.ref_viscosity/E->monitor.time_scale); /* scaling stress */
+
+ if(E->parallel.me == 0)fprintf(stderr, "ok22\n");
(E->problem_settings) (E);
-
+ if(E->parallel.me == 0)fprintf(stderr, "ok23\n");
viscosity_parameters(E);
return;
Modified: mc/3D/CitcomCU/trunk/src/Makefile.gzdir
===================================================================
--- mc/3D/CitcomCU/trunk/src/Makefile.gzdir 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Makefile.gzdir 2010-09-16 05:57:34 UTC (rev 17197)
@@ -18,15 +18,15 @@
#COMPRESS=/usr/bin/compress
COMPRESS=/bin/gzip
-GMT_DIR = ../../../GMTdev/GMT4.4.0/x86_64/lib/
+GMT_DIR = ../../../GMT4.5.3/
NETCDF_DIR = ../../../netcdf-3.6.2/
HC_DIR = ../../../hc-svn/
-LIB_PATH= -L$(HC_DIR)/objects/x86_64/ -L$(GMT_DIR) -L$(NETCDF_DIR)/lib/
+LIB_PATH= -L$(HC_DIR)/objects/x86_64/ -L$(GMT_DIR)/lib/ -L$(NETCDF_DIR)/lib/
#LIB_LIST= -lmpi
LIB_LIST= -lz -lggrd -lhc -lgmt -lnetcdf
LIB= $(LIB_PATH) $(LIB_LIST) -lm
-DEFINES = -DUSE_GZDIR -DUSE_GGRD -DUSE_GMT4 -I$(HC_DIR)/\
+DEFINES = -DUSE_GZDIR -DUSE_GGRD -I$(HC_DIR)/ -Werror \
-I$(GMT_DIR)/include/ -I$(NETCDF_DIR)/include/\
-DPREM_MODEL_FILE=\"$(PWD)/prem/prem.dat\"
Added: mc/3D/CitcomCU/trunk/src/Makefile.gzdir_ani
===================================================================
--- mc/3D/CitcomCU/trunk/src/Makefile.gzdir_ani (rev 0)
+++ mc/3D/CitcomCU/trunk/src/Makefile.gzdir_ani 2010-09-16 05:57:34 UTC (rev 17197)
@@ -0,0 +1,321 @@
+#
+# Makefile for "CitcomCU" allowing considerable variations
+#
+
+#
+# NOTES:
+# 1) Macro substitution is different for different version of make, for SGI,
+# CRAY, HP, etc. It seems necessary to go through and replace the
+# CEXT ... FOBJFLAG macros manually.
+#
+# 2) Put FLAGS=$(YOUR_MACHINE_FLAGS) chosen from the list that follows, add
+# machines, where necessary, try and find a better way than just defining
+# "const" to be null.
+#
+##############################################################################
+# Individual Machine Variations
+##############################################################################
+
+#COMPRESS=/usr/bin/compress
+COMPRESS=/bin/gzip
+GMT_DIR = ../../../GMT4.5.3/
+NETCDF_DIR = ../../../netcdf-3.6.2/
+HC_DIR = ../../../hc-svn/
+LIB_PATH= -L$(HC_DIR)/objects/x86_64/ -L$(GMT_DIR)/lib/ -L$(NETCDF_DIR)/lib/
+#LIB_LIST= -lmpi
+LIB_LIST= -lz -lggrd -lhc -lgmt -lnetcdf
+LIB= $(LIB_PATH) $(LIB_LIST) -lm
+
+DEFINES = -DUSE_GZDIR -DUSE_GGRD -DCITCOM_ALLOW_ANISOTROPIC_VISC \
+ -I$(HC_DIR)/ -Werror \
+ -I$(GMT_DIR)/include/ -I$(NETCDF_DIR)/include/\
+ -DPREM_MODEL_FILE=\"$(PWD)/prem/prem.dat\"
+
+###################################################################
+# Operating System Variables
+###################################################################
+
+#CC=/usr/local/mpich/bin/mpicc
+#CC=/opt/mpich/bin/mpicc
+#CC=/usr/lib/cmplrs/cc/gemc_cc
+#CC=$(HOME)/usr/local/mpich/bin/mpicc
+CC=mpicc
+F77=f77
+CPP=
+
+CEXT=c
+FEXT=F # which implies further action of cpp
+OBJEXT=o
+FOBJEXT=o
+OBJFLAG=-c
+FOBJFLAG=-c
+
+
+###################################
+# Choose your machine from here.
+###################################
+
+####################################
+# Dec Alpha, OSF 1
+AXPFLAGS= -unsigned -non_shared \
+ -math_library=fast -prefix=all -reentrancy=none -assume=noaccuracy_sensitive \
+ -unsigned_char -extern=strict_refdef -trapuv # -D_INTRINSICS
+AXPLDFLAGS= -unsigned -assume=noaccuracy_sensitive -non_shared -D_INTRINSICS
+AXPOPTIM= -O -Ublas -float -Olimit 1000 # -cord -feedback citcom.feedback
+####################################
+
+####################################
+# CRAY Unicos systems
+CRAYFLAGS=
+CRAYLDFLAGS=
+CRAYOPTIM=
+####################################
+
+####################################
+#IBM AIX systems
+AIXFLAGS= -D__aix__
+AIXLDFLAGS=
+AIXOPTIM= -O2 -qarch=pwr -qopt=pwr -Ublas
+####################################
+
+####################################
+# SUNOS systems
+SUNFLAGS= -D__sunos__ -Dconst=""
+SUNLDFLAGS=
+SUNOPTIM=-O -fsingle
+####################################
+
+####################################
+# Solaris systems
+SOLARISFLAGS= -D__solaris -Dconst="" -I/opt/mpi/include
+SOLARISLDFLAGS=-fast -lsocket -lnsl -lthread
+SOLARISOPTIM=-fast -xO4 -dalign -xtarget=ultra -xarch=v8plus -fsingle
+####################################
+
+####################################
+#HP running HPUX
+HPUXFLAGS=-Dconst=""
+HPUXLDFLAGS=
+HPUXOPTIM=+O3 +Onolimit +Odataprefetch
+#HPUXOPTIM=+O4 +Onolimit +Odataprefetch +Ofastaccess
+#HPUXOPTIM3=+O3 +Onolimit +Odataprefetch
+####################################
+
+####################################
+# SGI with IRIX
+SGIFLAGS=
+SGILDFLAGS=
+SGIOPTIM=-O -fsingle
+####################################
+
+#LinuxFLAGS=-I/usr/local/mpi/include
+#LinuxFLAGS=-I/opt/mpi/include
+#LinuxFLAGS=-I$(HOME)/usr/include
+LinuxFLAGS=
+LinuxLDFLAGS=
+#LinuxOPTIM=-g
+LinuxOPTIM=-O2 $(DEFINES)
+
+####################################
+#PARAGON
+PARAGONFLAGS=-nx
+PARAGONLDFLAGS=
+PARAGONOPTIM=-O4
+####################################
+
+####################################################################
+# Choose the actual flags for your machine here ...
+####################################################################
+
+FLAGS= $(LinuxFLAGS) -DCOMPRESS_BINARY=\"$(COMPRESS)\"
+LDFLAGS= $(LinuxLDFLAGS)
+OPTIM= $(LinuxOPTIM)
+
+#FLAGS= $(SOLARISFLAGS) -DCOMPRESS_BINARY=\"$(COMPRESS)\"
+#LDFLAGS= $(SOLARISLDFLAGS)
+#OPTIM= $(SOLARISOPTIM)
+
+####################################################################
+# CitcomCU's files .....
+####################################################################
+
+CFILES= Advection_diffusion.c\
+ Anisotropic_viscosity.c\
+ Boundary_conditions.c\
+ Citcom.c\
+ Composition_adv.c\
+ Construct_arrays.c\
+ Convection.c\
+ Drive_solvers.c\
+ Element_calculations.c\
+ Geometry_cartesian.c\
+ General_matrix_functions.c\
+ Global_operations.c\
+ Stokes_flow_Incomp.c\
+ Instructions.c\
+ Nodal_mesh.c\
+ Output.c\
+ Pan_problem_misc_functions.c\
+ Parallel_related.c\
+ Parsing.c\
+ Phase_change.c\
+ Process_buoyancy.c\
+ Process_velocity.c\
+ Profiling.c\
+ Shape_functions.c\
+ Size_does_matter.c\
+ Solver_multigrid.c\
+ Solver_conj_grad.c\
+ Sphere_harmonics.c\
+ Topo_gravity.c\
+ Viscosity_structures.c\
+ Output_gzdir.c\
+ Ggrd_handling.c
+
+
+
+HEADER = element_definitions.h\
+ global_defs.h\
+ viscosity_descriptions.h\
+ anisotropic_viscosity.h\
+ advection.h
+
+
+FFILES=#Blas_lapack_interfaces.F
+
+OBJFILES=$(CFILES:.c=.o) $(FFILES:.f=.o)
+
+
+####################################################################
+# Makefile rules follow
+####################################################################
+
+default: citcom.mpi
+
+citcom.mpi: $(OBJFILES) $(HEADER) Makefile
+ $(CC) $(OPTIM) $(FLAGS) $(LDFLAGS) -o citcom.mpi \
+ $(OBJFILES) $(FFTLIB) $(LIB)
+
+clean:
+ rm -f *.o
+
+clean-all:
+ rm -f *.o citcom.mpi
+
+smaller:
+ compress $(CFILES)
+
+larger:
+ uncompress $(CFILES)
+
+
+####################################################################
+# File dependencies follow
+####################################################################
+
+global_defs.h: viscosity_descriptions.h advection.h Convection_variables.h
+
+# The following entries can probably be automated from $CFILES etc
+
+Advection_diffusion.o: $(HEADER) Advection_diffusion.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Advection_diffusion.c
+
+Anisotropic_viscosity.o: $(HEADER) Anisotropic_viscosity.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Anisotropic_viscosity.c
+#
+Boundary_conditions.o: $(HEADER) Boundary_conditions.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Boundary_conditions.c
+#
+Citcom.o: $(HEADER) Citcom.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Citcom.c
+#
+Composition_adv.o: $(HEADER) Composition_adv.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Composition_adv.c
+#
+Construct_arrays.o: $(HEADER) Construct_arrays.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Construct_arrays.c
+#
+Convection.o: $(HEADER) Convection.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Convection.c
+#
+Drive_solvers.o: $(HEADER) Drive_solvers.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Drive_solvers.c
+#
+Element_calculations.o: $(HEADER) Element_calculations.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Element_calculations.c
+
+General_matrix_functions.o: $(HEADER) General_matrix_functions.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) General_matrix_functions.c
+
+Geometry_cartesian.o: $(HEADER) Geometry_cartesian.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Geometry_cartesian.c
+
+Global_operations.o: $(HEADER) Global_operations.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Global_operations.c
+
+Instructions.o: $(HEADER) Instructions.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Instructions.c
+
+Nodal_mesh.o: $(HEADER) Nodal_mesh.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Nodal_mesh.c
+
+Output.o: $(HEADER) Output.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Output.c
+
+Pan_problem_misc_functions.o: $(HEADER) Pan_problem_misc_functions.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Pan_problem_misc_functions.c
+
+Parallel_related.o: $(HEADER) Parallel_related.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Parallel_related.c
+
+Parsing.o: $(HEADER) Parsing.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Parsing.c
+
+Phase_change.o: $(HEADER) Phase_change.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Phase_change.c
+
+Process_velocity.o: $(HEADER) Process_velocity.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Process_velocity.c
+
+Process_buoyancy.o: $(HEADER) Process_buoyancy.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Process_buoyancy.c
+
+Profiling.o: $(HEADER) Profiling.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Profiling.c
+
+Shape_functions.o: $(HEADER) Shape_functions.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Shape_functions.c
+
+Size_does_matter.o: $(HEADER) Size_does_matter.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Size_does_matter.c
+
+Solver_conj_grad.o: $(HEADER) Solver_conj_grad.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Solver_conj_grad.c
+
+Solver_multigrid.o: $(HEADER) Solver_multigrid.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Solver_multigrid.c
+
+Sphere_harmonics.o: $(HEADER) Sphere_harmonics.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Sphere_harmonics.c
+
+Stokes_flow_Incomp.o: $(HEADER) Stokes_flow_Incomp.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Stokes_flow_Incomp.c
+
+Topo_gravity.o: $(HEADER) Topo_gravity.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Topo_gravity.c
+
+Viscosity_structures.o: $(HEADER) Viscosity_structures.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Viscosity_structures.c
+
+Blas_lapack_interfaces.f: Blas_lapack_interfaces.F
+ $(CPP) $(OPTIM) -P Blas_lapack_interfaces.F Blas_lapack_interfaces.f
+ $(F77) $(OPTIM) $(FOPTIM) -c Blas_lapack_interfaces.f
+
+Output_gzdir.o: $(HEADER) Output_gzdir.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Output_gzdir.c
+
+
+Ggrd_handling.o: $(HEADER) Ggrd_handling.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Ggrd_handling.c
+
+
Modified: mc/3D/CitcomCU/trunk/src/Output_gzdir.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -56,6 +56,9 @@
#include "element_definitions.h"
#include "global_defs.h"
#include "prototypes.h"
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+#include "anisotropic_viscosity.h"
+#endif
static gzFile *safe_gzopen(char *,char *);
FILE *safe_fopen(char *,char *);
@@ -112,7 +115,7 @@
sprintf(output_file,"if [ ! -s %s/%d ];then mkdir -p %s/%d;fi 2> /dev/null",
E->control.data_file2,file_number,E->control.data_file2,file_number);
system(output_file);
- fprintf(stderr,"making directory: %s\n",output_file);
+ //fprintf(stderr,"making directory: %s\n",output_file);
}
/* and wait for the other jobs */
@@ -207,8 +210,7 @@
sprintf(output_file,"%s/vtk_ecor.%d.gz",E->control.data_file2,E->parallel.me);
gzout = safe_gzopen(output_file,"w");
if(E->control.Rsphere){ /* spherical */
- if(E->parallel.me == 0)
- fprintf(stderr, "converting spherical to cartesian vtk coords to %s \n",output_file);
+ //if(E->parallel.me == 0)fprintf(stderr, "converting spherical to cartesian vtk coords to %s \n",output_file);
for(i=1;i <= E->lmesh.nno;i++) {
rtp2xyz((float)E->SX[3][i],(float)E->SX[1][i],(float)E->SX[2][i],locx);
gzprintf(gzout,"%9.6f %9.6f %9.6f\n",locx[0],locx[1],locx[2]);
@@ -242,8 +244,7 @@
E->ien[i].node[7]+offset,E->ien[i].node[8]+offset);
}
gzclose(gzout);
- if(E->parallel.me == 0)
- fprintf(stderr, "writing element connectivity %s \n",output_file);
+ //if(E->parallel.me == 0)fprintf(stderr, "writing element connectivity %s \n",output_file);
/* end vtk branch */
}
/* end init loop */
@@ -361,6 +362,24 @@
gzprintf(gzout,"%.4e\n",E->VI[E->mesh.levmax][i]);
}
gzclose(gzout);
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity){
+ sprintf(output_file,"%s/%d/avisc.%d.%d.gz",
+ E->control.data_file2,file_number, E->parallel.me,file_number);
+ gzout=safe_gzopen(output_file,"w");
+ gzprintf(gzout,"%3d %7d\n",1,E->lmesh.nno);
+ for(i=1;i<=E->lmesh.nno;i++){
+ gzprintf(gzout,"%10.4e %10.4e %10.4e %10.4e\n",
+ E->VI2[E->mesh.levmax][i],
+ E->VIn1[E->mesh.levmax][i],
+ E->VIn2[E->mesh.levmax][i],
+ E->VIn3[E->mesh.levmax][i]);
+ }
+ gzclose(gzout);
+
+
+ }
+#endif
}
if(vtk_comp_out){
/*
@@ -538,7 +557,7 @@
}
}
-if(E->parallel.me==0)fprintf(stderr,"vel output done\n");
+ if(E->parallel.me==0)fprintf(stderr,"vel output done\n");
return;
}
Modified: mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -914,6 +914,18 @@
xout[1] = rst * sin((double)phi); /* y */
xout[2] = (double)r * cos((double)theta);
}
+
+
+void xyz2rtp(float x,float y,float z,float *rout)
+{
+ float tmp1,tmp2;
+ tmp1 = x*x + y*y;
+ tmp2 = tmp1 + z*z;
+ rout[0] = sqrt(tmp2); /* r */
+ rout[1] = atan2(sqrt(tmp1),z); /* theta */
+ rout[2] = atan2(y,x); /* phi */
+}
+
void myerror(char *message, struct All_variables *E)
{
E->control.verbose = 1;
Modified: mc/3D/CitcomCU/trunk/src/Parallel_related.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Parallel_related.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Parallel_related.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -1095,6 +1095,7 @@
for(i = 1; i <= E->mesh.nsd; i++)
{
+
idb = 0;
for(k = 1; k <= 2; k++)
if(E->parallel.NUM_PASS[levmax].bound[i][k])
Modified: mc/3D/CitcomCU/trunk/src/Solver_multigrid.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Solver_multigrid.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Solver_multigrid.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -408,34 +408,72 @@
viscU = (float *)malloc((1 + E->lmesh.NNO[E->mesh.levmax]) * sizeof(float));
viscD = (float *)malloc((1 + vpts * E->lmesh.NNO[E->mesh.levmax - 1]) * sizeof(float));
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC /* allow for anisotropy */
+ if(E->viscosity.allow_anisotropic_viscosity){
+ /* isotropic */
+ for(lv = E->mesh.levmax; lv > E->mesh.levmin; lv--){
+ sl_minus = lv - 1;
+ if(E->viscosity.smooth_cycles == 1){
+ visc_from_gint_to_nodes(E, E->EVI[lv], viscU, lv);project_scalar(E, lv, viscU, viscD);visc_from_nodes_to_gint(E, viscD, E->EVI[sl_minus], sl_minus);
+ visc_from_gint_to_nodes(E, E->EVI2[lv], viscU, lv);project_scalar(E, lv, viscU, viscD);visc_from_nodes_to_gint(E, viscD, E->EVI2[sl_minus], sl_minus);
+ visc_from_gint_to_nodes(E, E->EVIn1[lv], viscU, lv);project_scalar(E, lv, viscU, viscD);visc_from_nodes_to_gint(E, viscD, E->EVIn1[sl_minus], sl_minus);
+ visc_from_gint_to_nodes(E, E->EVIn2[lv], viscU, lv);project_scalar(E, lv, viscU, viscD);visc_from_nodes_to_gint(E, viscD, E->EVIn2[sl_minus], sl_minus);
+ visc_from_gint_to_nodes(E, E->EVIn3[lv], viscU, lv);project_scalar(E, lv, viscU, viscD);visc_from_nodes_to_gint(E, viscD, E->EVIn3[sl_minus], sl_minus);
+ }else if(E->viscosity.smooth_cycles == 2){
+ visc_from_gint_to_ele(E, E->EVI[lv], viscU, lv);inject_scalar_e(E, lv, viscU, E->EVI[sl_minus]);
+ visc_from_gint_to_ele(E, E->EVI2[lv], viscU, lv);inject_scalar_e(E, lv, viscU, E->EVI2[sl_minus]);
+ visc_from_gint_to_ele(E, E->EVIn1[lv], viscU, lv);inject_scalar_e(E, lv, viscU, E->EVIn1[sl_minus]);
+ visc_from_gint_to_ele(E, E->EVIn2[lv], viscU, lv);inject_scalar_e(E, lv, viscU, E->EVIn2[sl_minus]);
+ visc_from_gint_to_ele(E, E->EVIn3[lv], viscU, lv);inject_scalar_e(E, lv, viscU, E->EVIn3[sl_minus]);
+ }else if(E->viscosity.smooth_cycles == 3){
+ visc_from_gint_to_ele(E, E->EVI[lv], viscU, lv);project_scalar_e(E, lv, viscU, viscD);visc_from_ele_to_gint(E, viscD, E->EVI[sl_minus], sl_minus);
+ visc_from_gint_to_ele(E, E->EVI2[lv], viscU, lv);project_scalar_e(E, lv, viscU, viscD);visc_from_ele_to_gint(E, viscD, E->EVI2[sl_minus], sl_minus);
+ visc_from_gint_to_ele(E, E->EVIn1[lv], viscU, lv);project_scalar_e(E, lv, viscU, viscD);visc_from_ele_to_gint(E, viscD, E->EVIn1[sl_minus], sl_minus);
+ visc_from_gint_to_ele(E, E->EVIn2[lv], viscU, lv);project_scalar_e(E, lv, viscU, viscD);visc_from_ele_to_gint(E, viscD, E->EVIn2[sl_minus], sl_minus);
+ visc_from_gint_to_ele(E, E->EVIn3[lv], viscU, lv);project_scalar_e(E, lv, viscU, viscD);visc_from_ele_to_gint(E, viscD, E->EVIn3[sl_minus], sl_minus);
+ }else if(E->viscosity.smooth_cycles == 0){
+ visc_from_gint_to_nodes(E, E->EVI[lv], viscU, lv);inject_scalar(E, lv, viscU, viscD);visc_from_nodes_to_gint(E, viscD, E->EVI[sl_minus], sl_minus);
+ visc_from_gint_to_nodes(E, E->EVI2[lv], viscU, lv);inject_scalar(E, lv, viscU, viscD);visc_from_nodes_to_gint(E, viscD, E->EVI2[sl_minus], sl_minus);
+ visc_from_gint_to_nodes(E, E->EVIn1[lv], viscU, lv);inject_scalar(E, lv, viscU, viscD);visc_from_nodes_to_gint(E, viscD, E->EVIn1[sl_minus], sl_minus);
+ visc_from_gint_to_nodes(E, E->EVIn2[lv], viscU, lv);inject_scalar(E, lv, viscU, viscD);visc_from_nodes_to_gint(E, viscD, E->EVIn2[sl_minus], sl_minus);
+ visc_from_gint_to_nodes(E, E->EVIn3[lv], viscU, lv);inject_scalar(E, lv, viscU, viscD);visc_from_nodes_to_gint(E, viscD, E->EVIn3[sl_minus], sl_minus);
- for(lv = E->mesh.levmax; lv > E->mesh.levmin; lv--)
- {
- sl_minus = lv - 1;
- if(E->viscosity.smooth_cycles == 1)
+ }
+ normalize_director_at_gint(E,E->EVIn1[sl_minus],E->EVIn2[sl_minus],E->EVIn3[sl_minus],sl_minus);
+ }
+ }else{
+#endif
+ /* isotropic */
+ for(lv = E->mesh.levmax; lv > E->mesh.levmin; lv--)
+ {
+ sl_minus = lv - 1;
+ if(E->viscosity.smooth_cycles == 1)
{
- visc_from_gint_to_nodes(E, E->EVI[lv], viscU, lv);
- project_scalar(E, lv, viscU, viscD);
- visc_from_nodes_to_gint(E, viscD, E->EVI[sl_minus], sl_minus);
+ visc_from_gint_to_nodes(E, E->EVI[lv], viscU, lv);
+ project_scalar(E, lv, viscU, viscD);
+ visc_from_nodes_to_gint(E, viscD, E->EVI[sl_minus], sl_minus);
}
- else if(E->viscosity.smooth_cycles == 2)
+ else if(E->viscosity.smooth_cycles == 2)
{
- visc_from_gint_to_ele(E, E->EVI[lv], viscU, lv);
- inject_scalar_e(E, lv, viscU, E->EVI[sl_minus]);
+ visc_from_gint_to_ele(E, E->EVI[lv], viscU, lv);
+ inject_scalar_e(E, lv, viscU, E->EVI[sl_minus]);
}
- else if(E->viscosity.smooth_cycles == 3)
+ else if(E->viscosity.smooth_cycles == 3)
{
- visc_from_gint_to_ele(E, E->EVI[lv], viscU, lv);
- project_scalar_e(E, lv, viscU, viscD);
- visc_from_ele_to_gint(E, viscD, E->EVI[sl_minus], sl_minus);
+ visc_from_gint_to_ele(E, E->EVI[lv], viscU, lv);
+ project_scalar_e(E, lv, viscU, viscD);
+ visc_from_ele_to_gint(E, viscD, E->EVI[sl_minus], sl_minus);
}
- else if(E->viscosity.smooth_cycles == 0)
+ else if(E->viscosity.smooth_cycles == 0)
{
- visc_from_gint_to_nodes(E, E->EVI[lv], viscU, lv);
- inject_scalar(E, lv, viscU, viscD);
- visc_from_nodes_to_gint(E, viscD, E->EVI[sl_minus], sl_minus);
+ visc_from_gint_to_nodes(E, E->EVI[lv], viscU, lv);
+ inject_scalar(E, lv, viscU, viscD);
+ visc_from_nodes_to_gint(E, viscD, E->EVI[sl_minus], sl_minus);
}
+ }
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
}
+#endif
free((void *)viscU);
free((void *)viscD);
Modified: mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -105,10 +105,9 @@
double alpha, delta, s2dotAhat, r0dotr0, r1dotz1;
//double residual, initial_residual, last_residual, res_magnitude, v_res;
double residual, initial_residual, res_magnitude, v_res;
-
+ char message[500];
double time0, time;
static double timea;
- //float dpressure, dvelocity, tole_comp;
float dpressure, dvelocity;
//const int dims = E->mesh.nsd;
@@ -165,7 +164,8 @@
assemble_div_u(E, V, r1, lev);
- residual = initial_residual = sqrt(global_pdot(E, r1, r1, lev) / gnpno);
+ E->monitor.pdotp = sqrt(global_pdot(E, r1, r1, lev) / gnpno);
+ residual = initial_residual = E->monitor.pdotp;
E->monitor.vdotv = sqrt(global_vdot(E, V, V, lev) / gneq);
@@ -182,15 +182,9 @@
res_magnitude = residual;
- if(E->control.print_convergence && E->parallel.me == 0)
- {
- fprintf(E->fp, "AhatP (%03d) after %g sec %g sec with div/v=%.3e for step %d\n",
- count, CPU_time0() - time0, CPU_time0() - timea, E->monitor.incompressibility, E->monitor.solution_cycles);
- /**/ fprintf(stderr, "AhatP (%03d) after %g sec %g sec with div/v=%.3e for step %d\n",
- count, CPU_time0() - time0, CPU_time0() - timea, E->monitor.incompressibility, E->monitor.solution_cycles);
- /**/
- }
+ generate_log_message(count,time0,timea,dvelocity, dpressure,E);
+
/* while( (count < *steps_max) && (E->monitor.incompressibility >= E->control.tole_comp || dvelocity >= imp) ) {
*/
while((count < *steps_max) && (dpressure >= imp || dvelocity >= imp))
@@ -237,19 +231,22 @@
V[j] -= alpha * u1[j];
assemble_div_u(E, V, Ah, lev);
- E->monitor.vdotv = global_vdot(E, V, V, E->mesh.levmax);
+ /* this is how it was computed before */
+ E->monitor.vdotv = global_vdot(E, V, V, lev);
+ E->monitor.pdotp = global_pdot(E, P, P, lev);
+
E->monitor.incompressibility = sqrt((gneq / gnpno) * (1.0e-32 + global_pdot(E, Ah, Ah, lev) / (1.0e-32 + E->monitor.vdotv)));
- dpressure = alpha * sqrt(global_pdot(E, s2, s2, lev) / (1.0e-32 + global_pdot(E, P, P, lev)));
+ dpressure = alpha * sqrt(global_pdot(E, s2, s2, lev) / (1.0e-32 + E->monitor.pdotp));
dvelocity = alpha * sqrt(global_vdot(E, u1, u1, lev) / (1.0e-32 + E->monitor.vdotv));
+ /* keep the normalized versions for the message */
+ E->monitor.vdotv = sqrt(E->monitor.vdotv/gneq);
+ E->monitor.pdotp = sqrt(E->monitor.pdotp/gnpno);
count++;
- if(E->control.print_convergence && E->parallel.me == 0)
- {
- fprintf(E->fp, "AhatP (%03d) after %g sec with div/v=%.3e, dv/v=%.3e & dp/p=%.3e for step %d\n", count, CPU_time0() - time0, E->monitor.incompressibility, dvelocity, dpressure, E->monitor.solution_cycles);
- /**/ fprintf(stderr, "AhatP (%03d) after %g sec with div/v=%.3e, dv/v=%.3e & dp/p=%.3e for step %d\n", count, CPU_time0() - time0, E->monitor.incompressibility, dvelocity, dpressure, E->monitor.solution_cycles);
- /**/ fflush(E->fp);
- }
+ generate_log_message(count,time0,timea,dvelocity, dpressure,E);
+
+
shuffle = s1;
s1 = s2;
s2 = shuffle;
@@ -296,7 +293,33 @@
return (residual);
}
+void generate_log_message(int count,double time0,double timea,double dvelocity,double dpressure, struct All_variables *E){
+ const int old_version=0;
+ char message[500];
+ if(E->control.print_convergence && E->parallel.me == 0){
+ if(old_version){ /* old version for backward compat */
+ if(count == 0)
+ /* first message was like this */
+ sprintf(message, "AhatP (%03d) after %g sec %g sec with div/v=%.3e for step %d\n",
+ count, CPU_time0() - time0, CPU_time0() - timea, E->monitor.incompressibility, E->monitor.solution_cycles);
+ else
+ /* other messages like that */
+ sprintf(message, "AhatP (%03d) after %g sec with div/v=%.3e, dv/v=%.3e & dp/p=%.3e for step %d\n",
+ count, CPU_time0() - time0, E->monitor.incompressibility, dvelocity, dpressure, E->monitor.solution_cycles);
+ }else{
+ sprintf(message, "AhatP (%03d) after %8.2f sec with v %12.6e dv/v %9.3e div/v %9.3e p %12.6e dp/p %9.3e for step %d\n",
+ count, CPU_time0() - time0,
+ E->monitor.vdotv,E->monitor.incompressibility, dvelocity,
+ E->monitor.pdotp,dpressure, E->monitor.solution_cycles);
+
+
+ }
+ fprintf(E->fp,"%s",message);//fflush(E->fp);
+ fprintf(stderr,"%s",message);
+ }
+}
+
/* ========================================================================== */
void v_from_vector(struct All_variables *E, float **V, double *F)
Modified: mc/3D/CitcomCU/trunk/src/Topo_gravity.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Topo_gravity.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Topo_gravity.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -41,6 +41,11 @@
#include "element_definitions.h"
#include "global_defs.h"
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+#include "anisotropic_viscosity.h"
+#endif
+
+
#define c_re(a) a.real
#define c_im(a) a.imag
typedef struct compl
@@ -93,7 +98,12 @@
const int elz = E->lmesh.elz;
//const int ely = E->lmesh.ely;
const int lev = E->mesh.levmax;
-
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity){
+ fprintf(stderr,"anisotropic viscosity not implemented for CBF topo\n");
+ parallel_process_termination();
+ }
+#endif
lnsf = E->lmesh.nsf;
eltTU = (float *)malloc((1 + Tsize) * sizeof(float));
@@ -304,6 +314,10 @@
//float pre[9], el_volume, tww[9], Visc, a, b;
float pre[9];
double rtf[4][9];
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ double D[6][6],n[3],eps[6],str[6];
+ int l1,l2;
+#endif
const int dims = E->mesh.nsd;
//const int dofs = E->mesh.dof;
@@ -339,7 +353,8 @@
Sxy = 0.0;
Sxz = 0.0;
Szy = 0.0;
- get_rtf(E, e, 0, rtf, lev);
+ if(E->control.Rsphere)
+ get_rtf(E, e, 0, rtf, lev);
for(j = 1; j <= vpts; j++)
{
@@ -371,7 +386,7 @@
Vxy[i] += (VX[j] * E->gNX[e].vpt[GNVXINDEX(1, j, i)] + VY[j] * E->gNX[e].vpt[GNVXINDEX(0, j, i)]);
Vzy[i] += (VY[j] * E->gNX[e].vpt[GNVXINDEX(2, j, i)] + VZ[j] * E->gNX[e].vpt[GNVXINDEX(1, j, i)]);
}
- else if(E->control.Rsphere)
+ else if(E->control.Rsphere) /* inaccrurate */
for(j = 1; j <= ends; j++)
{
Vzz[i] += VZ[j] * E->gNX[e].vpt[GNVXINDEX(2, j, i)];
@@ -382,14 +397,40 @@
Vzy[i] += VY[j] * E->gNX[e].vpt[GNVXINDEX(2, j, i)] + rtf[3][i] * (VZ[j] * E->gNX[e].vpt[GNVXINDEX(1, j, i)] / sin(rtf[1][i]) - VY[j] * E->N.vpt[GNVINDEX(j, i)]);
}
-
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity){ /* general anisotropic */
+
+ l1 = (e-1)*vpts+i;
+ get_constitutive(D,E->mesh.levmax,l1,rtf[1][i],rtf[2][i],(E->control.Rsphere),E);
+ eps[0] = Vxx[i];
+ eps[1] = Vyy[i];
+ eps[2] = Vzz[i];
+ eps[3] = Vxy[i];
+ eps[4] = Vxz[i];
+ eps[5] = Vzy[i];
+ for(l1=0;l1 < 6;l1++){
+ str[l1]=0.0;
+ for(l2=0;l2 < 6;l2++)
+ str[l1] += D[l1][l2] * eps[l2];
+ }
+ Sxx += pre[i] * str[0];
+ Syy += pre[i] * str[1];
+ Szz += pre[i] * str[2];
+ Sxy += pre[i] * str[3];
+ Sxz += pre[i] * str[4];
+ Szy += pre[i] * str[5];
+ }else{
+#endif
Sxx += 2.0 * pre[i] * Vxx[i];
Syy += 2.0 * pre[i] * Vyy[i];
Szz += 2.0 * pre[i] * Vzz[i];
Sxy += pre[i] * Vxy[i];
Sxz += pre[i] * Vxz[i];
Szy += pre[i] * Vzy[i];
- }
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ }
+#endif
+ } /* end i-vtps loop */
Sxx /= E->eco[e].area;
Modified: mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2010-09-16 05:57:34 UTC (rev 17197)
@@ -45,6 +45,9 @@
#include <string.h>
#include "element_definitions.h"
#include "global_defs.h"
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+#include "anisotropic_viscosity.h"
+#endif
static void visc_from_B(struct All_variables *, float *, float *, int );
static void visc_from_C(struct All_variables *, float *, float *, int );
@@ -103,6 +106,49 @@
input_int("num_mat", &(E->viscosity.num_mat), "1",m);
+ /* four layers */
+ E->viscosity.zlm = 1.0;
+ E->viscosity.z410 = 1.0;
+ E->viscosity.zlith = 0.0;
+
+ E->viscosity.zbase_layer[0] = E->viscosity.zbase_layer[1] = -999;
+ if(E->control.CART3D) /* defaults could be betters */
+ {
+ input_float("z_lmantle", &(E->viscosity.zlm), "1.0", m);
+ input_float("z_410", &(E->viscosity.z410), "1.0", m);
+ input_float("z_lith", &(E->viscosity.zlith), "0.0", m);
+ input_float_vector("z_layer",E->viscosity.num_mat,(E->viscosity.zbase_layer),m);
+ }
+ else if(E->control.Rsphere)
+ {
+ input_float("r_lmantle", &(E->viscosity.zlm), "1.0", m);
+ input_float("r_410", &(E->viscosity.z410), "1.0", m);
+ input_float("r_lith", &(E->viscosity.zlith), "0.0", m);
+ input_float_vector("r_layer",E->viscosity.num_mat,(E->viscosity.zbase_layer),m);
+ }
+
+ /* no z_layer input found */
+ if((fabs(E->viscosity.zbase_layer[0]+999) < 1e-5) &&
+ (fabs(E->viscosity.zbase_layer[1]+999) < 1e-5)) {
+
+ if(E->viscosity.num_mat != 4)
+ myerror("error: either use z_layer for non dim layer depths, or set num_mat to four",E);
+
+ E->viscosity.zbase_layer[0] = E->viscosity.zlith;
+ E->viscosity.zbase_layer[1] = E->viscosity.z410;
+ E->viscosity.zbase_layer[2] = E->viscosity.zlm;
+ E->viscosity.zbase_layer[3] = 0.55;
+ }
+
+
+
+
+
+
+
+
+
+
input_float_vector("viscT", E->viscosity.num_mat, (E->viscosity.T),m); /* redundant */
input_float_vector("viscT1", E->viscosity.num_mat, (E->viscosity.T),m);
input_float_vector("viscZ", E->viscosity.num_mat, (E->viscosity.Z),m);
@@ -137,7 +183,38 @@
/*
*/
-
+ 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_anisotropic_viscosity){ /* read additional
+ parameters for
+ anisotropic
+ viscosity */
+ 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
+ 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
+ */
+
+ input_boolean("anivisc_start_from_iso",
+ &(E->viscosity.anivisc_start_from_iso),"on",m); /* start
+ from
+ isotropic
+ solution? */
+ }
+#endif
+
+
+
+
/* composition factors */
input_float_vector("pre_comp",2,(E->viscosity.pre_comp),m);
@@ -215,8 +292,16 @@
const int vpts = vpoints[E->mesh.nsd];
-
+#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);
+ }
+#endif
+
if(E->viscosity.TDEPV)
visc_from_T(E, visc, evisc, propogate);
else
@@ -233,8 +318,9 @@
if(E->viscosity.SMOOTH)
- apply_viscosity_smoother(E, visc, evisc);
+ apply_viscosity_smoother(E, visc, evisc);
+
if(E->viscosity.MAX)
{
for(i = 1; i <= E->lmesh.nel; i++)
@@ -256,7 +342,16 @@
#ifdef USE_GZDIR
/* this is much preferred over v_to_nodes */
visc_from_gint_to_nodes(E,evisc,visc,E->mesh.levmax);
+#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);
+ visc_from_gint_to_nodes(E,E->EVIn3[E->mesh.levmax], E->VIn3[E->mesh.levmax],E->mesh.levmax);
+ normalize_director_at_nodes(E,E->VIn1[E->mesh.levmax],E->VIn2[E->mesh.levmax],E->VIn3[E->mesh.levmax],E->mesh.levmax);
+ }
#endif
+#endif
/* v_to_nodes(E,evisc,visc,E->mesh.levmax); */
return;
@@ -275,6 +370,18 @@
{
p_to_centres(E, visc, ViscCentre, E->mesh.levmax);
p_to_nodes(E, ViscCentre, visc, E->mesh.levmax);
+
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC /* allow for anisotropy */
+ if(E->viscosity.allow_anisotropic_viscosity){
+ p_to_centres(E, E->EVI2[E->mesh.levmax], ViscCentre, E->mesh.levmax);p_to_nodes(E, ViscCentre, E->EVI2[E->mesh.levmax], E->mesh.levmax);
+ p_to_centres(E, E->EVIn1[E->mesh.levmax], ViscCentre, E->mesh.levmax);p_to_nodes(E, ViscCentre, E->EVIn1[E->mesh.levmax], E->mesh.levmax);
+ p_to_centres(E, E->EVIn2[E->mesh.levmax], ViscCentre, E->mesh.levmax);p_to_nodes(E, ViscCentre, E->EVIn2[E->mesh.levmax], E->mesh.levmax);
+ p_to_centres(E, E->EVIn3[E->mesh.levmax], ViscCentre, E->mesh.levmax);p_to_nodes(E, ViscCentre, E->EVIn3[E->mesh.levmax], E->mesh.levmax);
+ normalize_director_at_gint(E,E->EVIn1[E->mesh.levmax],E->EVIn2[E->mesh.levmax],E->EVIn3[E->mesh.levmax],E->mesh.levmax);
+ }
+#endif
+
+
}
free((void *)ViscCentre);
@@ -706,23 +813,16 @@
int layers(struct All_variables *E, float x3)
{
int llayers = 0;
-
- /* this is the old logic, llayer = 4 would never get assigned */
- /* if(x3 >= E->viscosity.zlith) */
- /* llayers = 1; */
- /* else if((x3 < E->viscosity.zlith) && (x3 >= E->viscosity.zlm)) */
- /* llayers = 2; */
- /* else if(x3 < E->viscosity.zlm) */
- /* llayers = 3; */
-
- if(x3 >= E->viscosity.zlith) /* above 410 */
- llayers = 1;
- else if(x3 >= E->viscosity.z410) /* above 410 */
- llayers = 2;
- else if(x3 >= E->viscosity.zlm) /* above 660 */
- llayers = 3;
- else /* lower mantle */
- llayers = 4;
+ int i,ncheck;
+ ncheck = E->viscosity.num_mat-1;
+ for(i=0;i < ncheck;i++)
+ if(x3 >= E->viscosity.zbase_layer[i]){
+ llayers = i+1;
+ break;
+ }
+ if(!llayers)
+ llayers = E->viscosity.num_mat;
+
return (llayers);
}
Added: mc/3D/CitcomCU/trunk/src/anisotropic_viscosity.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/anisotropic_viscosity.h (rev 0)
+++ mc/3D/CitcomCU/trunk/src/anisotropic_viscosity.h 2010-09-16 05:57:34 UTC (rev 17197)
@@ -0,0 +1,34 @@
+#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 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 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 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 , double , double ,
+ int,struct All_variables *);
+
+#define __CITCOM_READ_ANIVISC_HEADER__
+#endif
Modified: mc/3D/CitcomCU/trunk/src/global_defs.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/global_defs.h 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/global_defs.h 2010-09-16 05:57:34 UTC (rev 17197)
@@ -648,6 +648,7 @@
int sobs_output_cols;
int solution_cycles;
+ int visc_iter_count;
float time_scale,time_scale_ma;
float length_scale;
@@ -659,13 +660,14 @@
float tpgscale;
float grvscale;
+
float delta_v_last_soln;
float elapsed_time;
float elapsed_time_vsoln;
float elapsed_time_vsoln1;
float reference_stress;
float incompressibility;
- float vdotv;
+ float vdotv,pdotp;
float nond_av_heat_fl;
float nond_av_adv_hfl;
float cpu_time_elapsed;
@@ -748,6 +750,7 @@
int stokes;
+
int check_t_irange,check_c_irange;
float Ra_670, clapeyron670, transT670, width670;
@@ -976,6 +979,15 @@
float *VB[4], *TB[4]; /* boundary conditions for V,T defined everywhere */
float *TW[MAX_LEVELS]; /* nodal weightings */
+
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ float *VI2[MAX_LEVELS],*EVI2[MAX_LEVELS];
+ float *VIn1[MAX_LEVELS],*EVIn1[MAX_LEVELS];
+ float *VIn2[MAX_LEVELS],*EVIn2[MAX_LEVELS];
+ float *VIn3[MAX_LEVELS],*EVIn3[MAX_LEVELS];
+#endif
+
+
int *surf_node, *surf_element;
int *mat; /* properties of mat */
unsigned int *node; /* properties of node */
Modified: mc/3D/CitcomCU/trunk/src/prototypes.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/prototypes.h 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/prototypes.h 2010-09-16 05:57:34 UTC (rev 17197)
@@ -1,319 +1,318 @@
/* Advection_diffusion.c */
-void advection_diffusion_parameters(struct All_variables *E);
-void advection_diffusion_allocate_memory(struct All_variables *E);
-void PG_timestep_particle(struct All_variables *E);
-void PG_timestep(struct All_variables *E);
-void predictor(struct All_variables *E, float *field, float *fielddot);
-void corrector(struct All_variables *E, float *field, float *fielddot, float *Dfielddot);
-void pg_solver(struct All_variables *E, float *T, float *Tdot, float *DTdot, float **V, struct SOURCES Q0, float diff, int bc, float **TBC, unsigned int *FLAGS);
-void pg_shape_fn(struct All_variables *E, int el, struct Shape_function *PG, float **V, double rtf[4][9], float diffusion);
-void element_residual(struct All_variables *E, int el, struct Shape_function PG, float **vel, float *field, float *fielddot, struct SOURCES Q0, double Eres[9], double rtf[4][9], float diff, float **BC, unsigned int *FLAGS);
-void std_timestep(struct All_variables *E);
-void process_heating(struct All_variables *E);
+void advection_diffusion_parameters(struct All_variables *);
+void advection_diffusion_allocate_memory(struct All_variables *);
+void PG_timestep_particle(struct All_variables *);
+void PG_timestep(struct All_variables *);
+void predictor(struct All_variables *, float *, float *);
+void corrector(struct All_variables *, float *, float *, float *);
+void pg_solver(struct All_variables *, float *, float *, float *, float **, struct SOURCES, float, int, float **, unsigned int *);
+void pg_shape_fn(struct All_variables *, int, struct Shape_function *, float **, double [4][9], float);
+void element_residual(struct All_variables *, int, struct Shape_function, float **, float *, float *, struct SOURCES, double [9], double [4][9], float, float **, unsigned int *);
+void std_timestep(struct All_variables *);
+void process_heating(struct All_variables *);
+/* Anisotropic_viscosity.c */
/* Boundary_conditions.c */
-void velocity_boundary_conditions(struct All_variables *E);
-void temperature_boundary_conditions(struct All_variables *E);
-void velocity_refl_vert_bc(struct All_variables *E);
-void temperature_refl_vert_bc(struct All_variables *E);
-void temperature_imposed_botm_bcs(struct All_variables *E, float *BC[], int dirn);
-void horizontal_bc(struct All_variables *E, float *BC[], int ROW, int dirn, float value, unsigned int mask, char onoff, int level);
-void velocity_apply_periodic_bcs(struct All_variables *E);
-void temperature_apply_periodic_bcs(struct All_variables *E);
-void strip_bcs_from_residual(struct All_variables *E, double *Res, int level);
-void temperatures_conform_bcs(struct All_variables *E);
-void velocities_conform_bcs(struct All_variables *E, double *U);
-void equalize_id_ien_lm(struct All_variables *E);
+void velocity_boundary_conditions(struct All_variables *);
+void temperature_boundary_conditions(struct All_variables *);
+void velocity_refl_vert_bc(struct All_variables *);
+void temperature_refl_vert_bc(struct All_variables *);
+void temperature_imposed_botm_bcs(struct All_variables *, float *[], int);
+void horizontal_bc(struct All_variables *, float *[], int, int, float, unsigned int, char, int);
+void velocity_apply_periodic_bcs(struct All_variables *);
+void temperature_apply_periodic_bcs(struct All_variables *);
+void strip_bcs_from_residual(struct All_variables *, double *, int);
+void temperatures_conform_bcs(struct All_variables *);
+void velocities_conform_bcs(struct All_variables *, double *);
+void equalize_id_ien_lm(struct All_variables *);
/* Citcom.c */
-int main(int argc, char **argv);
+int main(int, char **);
/* Composition_adv.c */
-void Runge_Kutta(struct All_variables *E, float *C, float *V[4], int on_off);
-void Euler(struct All_variables *E, float *C, float *V[4], int on_off);
-void transfer_markers_processors(struct All_variables *E, int on_off);
-void unify_markers_array(struct All_variables *E, int no_tran, int no_recv);
-void prepare_transfer_arrays(struct All_variables *E);
-int locate_processor(struct All_variables *E, double XMC1, double XMC2, double XMC3);
-void get_C_from_markers(struct All_variables *E, float *C);
-void element_markers(struct All_variables *E, int con);
-void velocity_markers(struct All_variables *E, float *V[4], int con);
-int get_element(struct All_variables *E, double XMC1, double XMC2, double XMC3, double dX[4]);
-int in_the_domain(struct All_variables *E, double r, double t, double f);
-float area_of_4node1(float x1, float y1, float x2, float y2, float x3, float y3, float x4, float y4);
-float area_of_3node(float x1, float y1, float x2, float y2, float x3, float y3);
-float mean_of_3node(int a, float x1, float y1, float x2, float y2, float x3, float y3);
-float mean_of_4node(int a, float x1, float y1, float x2, float y2, float x3, float y3, float x4, float y4);
-float mean_of_5node(int a, float x1, float y1, float x2, float y2, float x3, float y3, float x4, float y4, float x5, float y5);
-float dist1(float XO[4], float XN[4]);
+void Runge_Kutta(struct All_variables *, float *, float *[4], int);
+void Euler(struct All_variables *, float *, float *[4], int);
+void transfer_markers_processors(struct All_variables *, int);
+void unify_markers_array(struct All_variables *, int, int);
+void prepare_transfer_arrays(struct All_variables *);
+int locate_processor(struct All_variables *, double, double, double);
+void get_C_from_markers(struct All_variables *, float *);
+void element_markers(struct All_variables *, int);
+void velocity_markers(struct All_variables *, float *[4], int);
+int get_element(struct All_variables *, double, double, double, double [4]);
+int in_the_domain(struct All_variables *, double, double, double);
+float area_of_4node1(float, float, float, float, float, float, float, float);
+float area_of_3node(float, float, float, float, float, float);
+float mean_of_3node(int, float, float, float, float, float, float);
+float mean_of_4node(int, float, float, float, float, float, float, float, float);
+float mean_of_5node(int, float, float, float, float, float, float, float, float, float, float);
+float dist1(float [4], float [4]);
/* Construct_arrays.c */
-void construct_ien(struct All_variables *E);
-void construct_id(struct All_variables *E);
-void construct_lm(struct All_variables *E);
-void construct_node_maps(struct All_variables *E);
-void construct_node_ks(struct All_variables *E);
-void construct_masks(struct All_variables *E);
-void construct_sub_element(struct All_variables *E);
-void construct_elt_ks(struct All_variables *E);
-void construct_elt_gs(struct All_variables *E);
-void construct_mat_group(struct All_variables *E);
-void construct_stiffness_B_matrix(struct All_variables *E);
-void rebuild_BI_on_boundary(struct All_variables *E);
+void construct_ien(struct All_variables *);
+void construct_id(struct All_variables *);
+void construct_lm(struct All_variables *);
+void construct_node_maps(struct All_variables *);
+void construct_node_ks(struct All_variables *);
+void construct_masks(struct All_variables *);
+void construct_sub_element(struct All_variables *);
+void construct_elt_ks(struct All_variables *);
+void construct_elt_gs(struct All_variables *);
+void construct_mat_group(struct All_variables *);
+void construct_stiffness_B_matrix(struct All_variables *);
+void rebuild_BI_on_boundary(struct All_variables *);
/* Convection.c */
-void set_convection_defaults(struct All_variables *E);
-void read_convection_settings(struct All_variables *E);
-void convection_derived_values(struct All_variables *E);
-void convection_allocate_memory(struct All_variables *E);
-void convection_initial_fields(struct All_variables *E);
-void convection_boundary_conditions(struct All_variables *E);
-void convection_initial_temperature(struct All_variables *E);
-void process_restart_tc(struct All_variables *E);
-void convection_initial_markers1(struct All_variables *E);
-void convection_initial_markers(struct All_variables *E, int use_element_nodes_for_init_c);
-void setup_plume_problem(struct All_variables *E);
-void PG_process(struct All_variables *E, int ii);
+void set_convection_defaults(struct All_variables *);
+void read_convection_settings(struct All_variables *);
+void convection_derived_values(struct All_variables *);
+void convection_allocate_memory(struct All_variables *);
+void convection_initial_fields(struct All_variables *);
+void convection_boundary_conditions(struct All_variables *);
+void convection_initial_temperature(struct All_variables *);
+void process_restart_tc(struct All_variables *);
+void convection_initial_markers1(struct All_variables *);
+void convection_initial_markers(struct All_variables *, int);
+void setup_plume_problem(struct All_variables *);
+void PG_process(struct All_variables *, int);
/* Drive_solvers.c */
-void general_stokes_solver(struct All_variables *E);
+void general_stokes_solver(struct All_variables *);
+int need_to_iterate(struct All_variables *);
/* Element_calculations.c */
-void assemble_forces(struct All_variables *E, int penalty);
-void get_elt_k(struct All_variables *E, int el, double elt_k[24 * 24], int lev, int iconv);
-void assemble_del2_u(struct All_variables *E, double *u, double *Au, int level, int strip_bcs);
-void e_assemble_del2_u(struct All_variables *E, double *u, double *Au, int level, int strip_bcs);
-void n_assemble_del2_u(struct All_variables *E, double *u, double *Au, int level, int strip_bcs);
-void build_diagonal_of_K(struct All_variables *E, int el, double elt_k[24 * 24], int level);
-void build_diagonal_of_Ahat(struct All_variables *E);
-void assemble_div_u(struct All_variables *E, double *U, double *divU, int level);
-void assemble_grad_p(struct All_variables *E, double *P, double *gradP, int lev);
-double assemble_dAhatp_entry(struct All_variables *E, int e, int level);
-void get_elt_g(struct All_variables *E, int el, higher_precision elt_del[24][1], int lev);
-void get_elt_h(struct All_variables *E, int el, double elt_h[1], int penalty);
-void get_elt_f(struct All_variables *E, int el, double elt_f[24], int penalty, int bcs);
-void get_aug_k(struct All_variables *E, int el, double elt_k[24 * 24], int level);
+void assemble_forces(struct All_variables *, int);
+void get_elt_k(struct All_variables *, int, double [24 * 24], int, int);
+void assemble_del2_u(struct All_variables *, double *, double *, int, int);
+void e_assemble_del2_u(struct All_variables *, double *, double *, int, int);
+void n_assemble_del2_u(struct All_variables *, double *, double *, int, int);
+void build_diagonal_of_K(struct All_variables *, int, double [24 * 24], int);
+void build_diagonal_of_Ahat(struct All_variables *);
+void assemble_div_u(struct All_variables *, double *, double *, int);
+void assemble_grad_p(struct All_variables *, double *, double *, int);
+double assemble_dAhatp_entry(struct All_variables *, int, int);
+void get_elt_g(struct All_variables *, int, higher_precision [24][1], int);
+void get_elt_h(struct All_variables *, int, double [1], int);
+void get_elt_f(struct All_variables *, int, double [24], int, int);
+void get_aug_k(struct All_variables *, int, double [24 * 24], int);
/* General_matrix_functions.c */
-double **dmatrix(int nrl, int nrh, int ncl, int nch);
-float **fmatrix(int nrl, int nrh, int ncl, int nch);
-void dfree_matrix(double **m, int nrl, int nrh, int ncl, int nch);
-void ffree_matrix(float **m, int nrl, int nrh, int ncl, int nch);
-double *dvector(int nl, int nh);
-float *fvector(int nl, int nh);
-void dfree_vector(double *v, int nl, int nh);
-void ffree_vector(float *v, int nl, int nh);
-int *sivector(int nl, int nh);
-void sifree_vector(int *v, int nl, int nh);
-double pdot(struct All_variables *E, double *A, double *B, int lev);
-double pselfdot(struct All_variables *E, double *A);
-double vdot(struct All_variables *E, double *A, double *B, int level);
-double vselfdot(struct All_variables *E, double *A, int level);
-double vfselfdot(struct All_variables *E, float *A, int level);
-float fdot(float *A, float *B, int n1, int n2);
-float fselfdot(float *A, int n1, int n2);
-float dot(struct All_variables *E, float *A, float *B);
-float selfdot(struct All_variables *E, float *A);
-void dvcopy(double *A, double *B, int a, int b);
-void vcopy(float *A, float *B, int a, int b);
-void vprod(double *R, double *A, double *B, int a, int b);
-float fnmax(struct All_variables *E, float *A, int a, int b);
-int solve_del2_u(struct All_variables *E, double *d0, double *F, double acc, int high_lev, int ic);
-double multi_grid(struct All_variables *E, double *d1, double *F, double *Au, double acc, int hl);
-double conj_grad(struct All_variables *E, double *d0, double *F, double *Au, double acc, int *cycles, int level);
-void jacobi(struct All_variables *E, double *d0, double *F, double *Ad, double acc, int *cycles, int level, int guess);
-void element_gauss_seidel(struct All_variables *E, double *d0, double *F, double *Ad, double acc, int *cycles, int level, int guess);
-void gauss_seidel1(struct All_variables *E, double *d0, double *F, double *Ad, double acc, int *cycles, int level, int guess);
-void gauss_seidel(struct All_variables *E, double *d0, double *F, double *Ad, double acc, int *cycles, int level, int guess);
-void print_elt_k(struct All_variables *E, double a[24 * 24]);
-double cofactor(double A[4][4], int i, int j, int n);
-double determinant(double A[4][4], int n);
-float area_of_4node(float x1, float y1, float x2, float y2, float x3, float y3, float x4, float y4);
-double modified_plgndr_a(int l, int m, double t);
-double sqrt_multis(int jj, int ii);
-double multis(int ii);
-int int_multis(int ii);
-double plgndr_a(int l, int m, double t);
-double sphere_h(int l, int m, double t, double f, int ic);
+double **dmatrix(int, int, int, int);
+float **fmatrix(int, int, int, int);
+void dfree_matrix(double **, int, int, int, int);
+void ffree_matrix(float **, int, int, int, int);
+double *dvector(int, int);
+float *fvector(int, int);
+void dfree_vector(double *, int, int);
+void ffree_vector(float *, int, int);
+int *sivector(int, int);
+void sifree_vector(int *, int, int);
+double pdot(struct All_variables *, double *, double *, int);
+double pselfdot(struct All_variables *, double *);
+double vdot(struct All_variables *, double *, double *, int);
+double vselfdot(struct All_variables *, double *, int);
+double vfselfdot(struct All_variables *, float *, int);
+float fdot(float *, float *, int, int);
+float fselfdot(float *, int, int);
+float dot(struct All_variables *, float *, float *);
+float selfdot(struct All_variables *, float *);
+void dvcopy(double *, double *, int, int);
+void vcopy(float *, float *, int, int);
+void vprod(double *, double *, double *, int, int);
+float fnmax(struct All_variables *, float *, int, int);
+int solve_del2_u(struct All_variables *, double *, double *, double, int, int);
+double multi_grid(struct All_variables *, double *, double *, double *, double, int);
+double conj_grad(struct All_variables *, double *, double *, double *, double, int *, int);
+void jacobi(struct All_variables *, double *, double *, double *, double, int *, int, int);
+void element_gauss_seidel(struct All_variables *, double *, double *, double *, double, int *, int, int);
+void gauss_seidel1(struct All_variables *, double *, double *, double *, double, int *, int, int);
+void gauss_seidel(struct All_variables *, double *, double *, double *, double, int *, int, int);
+void print_elt_k(struct All_variables *, double [24 * 24]);
+double cofactor(double [4][4], int, int, int);
+double determinant(double [4][4], int);
+float area_of_4node(float, float, float, float, float, float, float, float);
+double modified_plgndr_a(int, int, double);
+double sqrt_multis(int, int);
+double multis(int);
+int int_multis(int);
+double plgndr_a(int, int, double);
+double sphere_h(int, int, double, double, int);
/* Geometry_cartesian.c */
-void set_2dc_defaults(struct All_variables *E);
-void set_2pt5dc_defaults(struct All_variables *E);
-void set_3ds_defaults(struct All_variables *E);
-void set_3dc_defaults(struct All_variables *E);
+void set_2pt5dc_defaults(struct All_variables *);
+void set_3ds_defaults(struct All_variables *);
+void set_3dc_defaults(struct All_variables *);
/* Ggrd_handling.c */
-void convection_initial_temperature_ggrd(struct All_variables *E);
/* Global_operations.c */
-void remove_horiz_ave(struct All_variables *E, float *X, float *H, int store_or_not);
-void return_horiz_sum(struct All_variables *E, float *X, float *H, int nn);
-void return_horiz_ave(struct All_variables *E, float *X, float *H);
-float return_bulk_value(struct All_variables *E, float *Z, float z_thld, int average);
-double global_vdot(struct All_variables *E, double *A, double *B, int lev);
-double global_pdot(struct All_variables *E, double *A, double *B, int lev);
-float global_tdot(struct All_variables *E, float *A, float *B, int lev);
-float Tmax(struct All_variables *E, float *T);
-float global_fmin(struct All_variables *E, float a);
-float global_fmax(struct All_variables *E, float a);
-void sum_across_depth_sph1(struct All_variables *E, float *sphc, float *sphs);
-void sum_across_surface(struct All_variables *E, float *data, int total);
-void sum_across_surf_sph1(struct All_variables *E, float *sphc, float *sphs);
-void gather_TG_to_me0(struct All_variables *E, float *TG);
-void propogator_down_process(struct All_variables *E, float *Tadi);
-double sum_across_depth(struct All_variables *E, double temp1);
+void return_horiz_sum(struct All_variables *, float *, float *, int);
+void return_horiz_ave(struct All_variables *, float *, float *);
+float return_bulk_value(struct All_variables *, float *, float, int);
+double global_vdot(struct All_variables *, double *, double *, int);
+double global_pdot(struct All_variables *, double *, double *, int);
+float global_tdot(struct All_variables *, float *, float *, int);
+float Tmax(struct All_variables *, float *);
+float global_fmin(struct All_variables *, float);
+float global_fmax(struct All_variables *, float);
+void sum_across_depth_sph1(struct All_variables *, float *, float *);
+void sum_across_surface(struct All_variables *, float *, int);
+void sum_across_surf_sph1(struct All_variables *, float *, float *);
+void gather_TG_to_me0(struct All_variables *, float *);
+void propogator_down_process(struct All_variables *, float *);
+double sum_across_depth(struct All_variables *, double);
/* Instructions.c */
-void read_instructions(struct All_variables *E, int argc, char **argv);
-void allocate_common_vars(struct All_variables *E);
-void interruption(int signal_number);
-void global_default_values(struct All_variables *E);
-void global_derived_values(struct All_variables *E);
-void read_initial_settings(struct All_variables *E);
-void check_bc_consistency(struct All_variables *E);
-void set_up_nonmg_aliases(struct All_variables *E);
-void report(struct All_variables *E, char *string);
-void record(struct All_variables *E, char *string);
-void common_initial_fields(struct All_variables *E);
-void initial_pressure(struct All_variables *E);
-void initial_velocity(struct All_variables *E);
+void read_instructions(struct All_variables *, int, char **);
+void allocate_common_vars(struct All_variables *);
+void interruption(int);
+void global_default_values(struct All_variables *);
+void global_derived_values(struct All_variables *);
+void read_initial_settings(struct All_variables *);
+void check_bc_consistency(struct All_variables *);
+void set_up_nonmg_aliases(struct All_variables *);
+void report(struct All_variables *, char *);
+void record(struct All_variables *, char *);
+void common_initial_fields(struct All_variables *);
+void initial_pressure(struct All_variables *);
+void initial_velocity(struct All_variables *);
/* Nodal_mesh.c */
-void node_locations(struct All_variables *E);
-void pre_interpolation(struct All_variables *E);
-void dlogical_mesh_to_real(struct All_variables *E, double *data, int level);
-void flogical_mesh_to_real(struct All_variables *E, float *data, int level);
-void p_to_nodes(struct All_variables *E, double *P, float *PN, int lev);
-void p_to_centres(struct All_variables *E, float *PN, double *P, int lev);
-void v_to_intpts(struct All_variables *E, float *VN, float *VE, int lev);
-void v_to_nodes(struct All_variables *E, float *VE, float *VN, int lev);
-void visc_to_intpts(struct All_variables *E, float *VN, float *VE, int lev);
-void visc_to_nodes(struct All_variables *E, float *VE, float *VN, int lev);
-void visc_from_ele_to_gint(struct All_variables *E, float *VN, float *VE, int lev);
-void visc_from_gint_to_ele(struct All_variables *E, float *VE, float *VN, int lev);
-void visc_from_gint_to_nodes(struct All_variables *E, float *VE, float *VN, int lev);
-void visc_from_nodes_to_gint(struct All_variables *E, float *VN, float *VE, int lev);
+void node_locations(struct All_variables *);
+void pre_interpolation(struct All_variables *);
+void dlogical_mesh_to_real(struct All_variables *, double *, int);
+void flogical_mesh_to_real(struct All_variables *, float *, int);
+void p_to_nodes(struct All_variables *, double *, float *, int);
+void p_to_centres(struct All_variables *, float *, double *, int);
+void v_to_intpts(struct All_variables *, float *, float *, int);
+void v_to_nodes(struct All_variables *, float *, float *, int);
+void visc_to_intpts(struct All_variables *, float *, float *, int);
+void visc_to_nodes(struct All_variables *, float *, float *, int);
+void visc_from_ele_to_gint(struct All_variables *, float *, float *, int);
+void visc_from_gint_to_ele(struct All_variables *, float *, float *, int);
+void visc_from_gint_to_nodes(struct All_variables *, float *, float *, int);
+void visc_from_nodes_to_gint(struct All_variables *, float *, float *, int);
/* Output.c */
-void output_velo_related(struct All_variables *E, int file_number);
-void output_velo_related_binary(struct All_variables *E, int file_number);
-void output_temp(struct All_variables *E, int file_number);
-void process_restart(struct All_variables *E);
-void print_field_spectral_regular(struct All_variables *E, float *TG, float *sphc, float *sphs, int proc_loc, char *filen);
+void output_velo_related_binary(struct All_variables *, int);
+void output_temp(struct All_variables *, int);
+void process_restart(struct All_variables *);
+void print_field_spectral_regular(struct All_variables *, float *, float *, float *, int, char *);
/* Output_gzdir.c */
-void output_velo_related_gzdir(struct All_variables *E, int file_number);
-void process_restart_tc_gzdir(struct All_variables *E);
+void output_velo_related_gzdir(struct All_variables *, int);
+void process_restart_tc_gzdir(struct All_variables *);
/* Pan_problem_misc_functions.c */
int get_process_identifier(void);
-void unique_copy_file(struct All_variables *E, char *name, char *comment);
-void thermal_buoyancy(struct All_variables *E);
-double SIN_D(double x);
-double COT_D(double x);
-void *Malloc1(int bytes, char *file, int line);
-int read_previous_field(struct All_variables *E, float *field, char *name, char *abbr);
-void fcopy_interpolating(struct All_variables *E, float *X, float *Z, float *Y, int nx, int nz, int ny, float *T, float *TT);
-float cross2d(float x11, float x12, float x21, float x22, int D);
-void field_arbitrary_rectangle_file(struct All_variables *E, int parse_and_apply, struct Rect *RECT, char *name, float *field, int BC, unsigned int *bcbitf, unsigned int bcmask_on, unsigned int bcmask_off);
-void field_arbitrary_rectangle(struct All_variables *E, struct Rect *RECT, float *field, int BC, unsigned int *bcbitf, unsigned int bcmask_on, unsigned int bcmask_off);
-void field_arbitrary_circle_file(struct All_variables *E, int parse_and_apply, struct Circ *CIRC, char *name, float *field, int BC, unsigned int *bcbitf, unsigned int bcmask_on, unsigned int bcmask_off);
-void field_arbitrary_circle(struct All_variables *E, struct Circ *CIRC, float *field, int BC, unsigned int *bcbitf, unsigned int bcmask_on, unsigned int bcmask_off);
-void field_arbitrary_harmonic_file(struct All_variables *E, int parse_and_apply, struct Harm *HARM, char *name, float *field, int BC, unsigned int *bcbitf, unsigned int bcmask_on, unsigned int bcmask_off);
-void field_arbitrary_harmonic(struct All_variables *E, struct Harm *HARM, float *field, int BC, unsigned int *bcbitf, unsigned int bcmask_on, unsigned int bcmask_off);
-double myatan(double y, double x);
-FILE *safe_fopen(char *name, char *mode);
-void *safe_malloc(size_t size);
-void calc_cbase_at_tp(float theta, float phi, float *base);
-void convert_pvec_to_cvec(float vr, float vt, float vp, float *base, float *cvec);
-void rtp2xyz(float r, float theta, float phi, float *xout);
-void myerror(char *message, struct All_variables *E);
+void unique_copy_file(struct All_variables *, char *, char *);
+void thermal_buoyancy(struct All_variables *);
+double SIN_D(double);
+double COT_D(double);
+void *Malloc1(int, char *, int);
+int read_previous_field(struct All_variables *, float *, char *, char *);
+void fcopy_interpolating(struct All_variables *, float *, float *, float *, int, int, int, float *, float *);
+float cross2d(float, float, float, float, int);
+void field_arbitrary_rectangle_file(struct All_variables *, int, struct Rect *, char *, float *, int, unsigned int *, unsigned int, unsigned int);
+void field_arbitrary_rectangle(struct All_variables *, struct Rect *, float *, int, unsigned int *, unsigned int, unsigned int);
+void field_arbitrary_circle_file(struct All_variables *, int, struct Circ *, char *, float *, int, unsigned int *, unsigned int, unsigned int);
+void field_arbitrary_circle(struct All_variables *, struct Circ *, float *, int, unsigned int *, unsigned int, unsigned int);
+void field_arbitrary_harmonic_file(struct All_variables *, int, struct Harm *, char *, float *, int, unsigned int *, unsigned int, unsigned int);
+void field_arbitrary_harmonic(struct All_variables *, struct Harm *, float *, int, unsigned int *, unsigned int, unsigned int);
+double myatan(double, double);
+FILE *safe_fopen(char *, char *);
+void *safe_malloc(size_t);
+void calc_cbase_at_tp(float, float, float *);
+void convert_pvec_to_cvec(float, float, float, float *, float *);
+void rtp2xyz(float, float, float, float *);
+void myerror(char *, struct All_variables *);
/* Parallel_related.c */
-void parallel_process_initilization(struct All_variables *E, int argc, char **argv);
void parallel_process_termination(void);
-void parallel_domain_decomp1(struct All_variables *E);
-void parallel_shuffle_ele_and_id(struct All_variables *E);
-void parallel_shuffle_ele_and_id_bc1(struct All_variables *E);
-void parallel_shuffle_ele_and_id_bc2(struct All_variables *E);
-void parallel_communication_routs(struct All_variables *E);
-void parallel_communication_routs1(struct All_variables *E);
-void parallel_communication_routs2(struct All_variables *E);
-void parallel_communication_routs3(struct All_variables *E);
-void parallel_communication_routs4(struct All_variables *E);
-void exchange_number_rec_markers(struct All_variables *E);
-void exchange_markers(struct All_variables *E);
-void exchange_id_d20(struct All_variables *E, double *U, int lev);
-void exchange_node_f20(struct All_variables *E, float *U, int lev);
+void parallel_domain_decomp1(struct All_variables *);
+void parallel_shuffle_ele_and_id(struct All_variables *);
+void parallel_shuffle_ele_and_id_bc1(struct All_variables *);
+void parallel_shuffle_ele_and_id_bc2(struct All_variables *);
+void parallel_communication_routs(struct All_variables *);
+void parallel_communication_routs1(struct All_variables *);
+void parallel_communication_routs2(struct All_variables *);
+void parallel_communication_routs3(struct All_variables *);
+void parallel_communication_routs4(struct All_variables *);
+void exchange_number_rec_markers(struct All_variables *);
+void exchange_markers(struct All_variables *);
+void exchange_id_d20(struct All_variables *, double *, int);
+void exchange_node_f20(struct All_variables *, float *, int);
double CPU_time0(void);
void parallel_process_sync(void);
/* Parsing.c */
-void setup_parser(struct All_variables *E, char *filename);
-void shutdown_parser(struct All_variables *E);
-void add_to_parameter_list(char *name, char *value);
-int compute_parameter_hash_table(char *s);
-int input_int(char *name, int *value, char *interpret, int m);
-int input_string(char *name, char *value, char *Default, int m);
-int input_boolean(char *name, int *value, char *interpret, int m);
-int input_float(char *name, float *value, char *interpret, int m);
-int input_double(char *name, double *value, char *interpret, int m);
-int input_int_vector(char *name, int number, int *value, int m);
-int input_char_vector(char *name, int number, char *value, int m);
-int input_float_vector(char *name, int number, float *value, int m);
-int input_double_vector(char *name, int number, double *value, int m);
-int interpret_control_string(char *interpret, int *essential, double *Default, double *minvalue, double *maxvalue);
+void setup_parser(struct All_variables *, char *);
+void shutdown_parser(struct All_variables *);
+void add_to_parameter_list(char *, char *);
+int compute_parameter_hash_table(char *);
+int input_int(char *, int *, char *, int);
+int input_string(char *, char *, char *, int);
+int input_boolean(char *, int *, char *, int);
+int input_float(char *, float *, char *, int);
+int input_double(char *, double *, char *, int);
+int input_int_vector(char *, int, int *, int);
+int input_char_vector(char *, int, char *, int);
+int input_float_vector(char *, int, float *, int);
+int input_double_vector(char *, int, double *, int);
+int interpret_control_string(char *, int *, double *, double *, double *);
/* Phase_change.c */
-void phase_change(struct All_variables *E, float *B6, float *B_b6, float *B4, float *B_b4);
/* Process_buoyancy.c */
-void process_temp_field(struct All_variables *E, int ii);
-void heat_flux(struct All_variables *E);
-void heat_flux1(struct All_variables *E);
-void plume_buoyancy_flux(struct All_variables *E);
+void heat_flux(struct All_variables *);
+void heat_flux1(struct All_variables *);
+void plume_buoyancy_flux(struct All_variables *);
/* Process_velocity.c */
-void process_new_velocity(struct All_variables *E, int ii);
-void get_surface_velo(struct All_variables *E, float *SV);
-void get_ele_visc(struct All_variables *E, float *EV);
-void get_surf_stress(struct All_variables *E, float *SXX, float *SYY, float *SZZ, float *SXY, float *SXZ, float *SZY);
-void averages(struct All_variables *E);
+void get_surface_velo(struct All_variables *, float *);
+void get_ele_visc(struct All_variables *, float *);
+void get_surf_stress(struct All_variables *, float *, float *, float *, float *, float *, float *);
+void averages(struct All_variables *);
/* Profiling.c */
float CPU_time(void);
/* Shape_functions.c */
-void construct_shape_functions(struct All_variables *E);
-double lpoly(int p, double y);
-double lpolydash(int p, double y);
+double lpoly(int, double);
+double lpolydash(int, double);
/* Size_does_matter.c */
-void twiddle_thumbs(struct All_variables *yawn, int scratch_groin);
-void get_global_shape_fn(struct All_variables *E, int el, int pressure, double rtf[4][9], int sphere, int level);
-void form_rtf_bc(int k, double x[4], double rtf[4][9], double bc[4][4]);
-void get_rtf(struct All_variables *E, int el, int pressure, double rtf[4][9], int lev);
-void construct_c3x3matrix_el(struct All_variables *E, int el, struct CC *cc, struct CCX *ccx, int lev, int pressure);
-void get_global_1d_shape_fn(struct All_variables *E, int el, struct Shape_function1 *GM, struct Shape_function1_dA *dGammax, int top);
-void get_global_1d_shape_fn1(struct All_variables *E, int el, struct Shape_function1 *GM, struct Shape_function1_dA *dGammax, int top);
-void mass_matrix(struct All_variables *E);
+void get_global_shape_fn(struct All_variables *, int, int, double [4][9], int, int);
+void form_rtf_bc(int, double [4], double [4][9], double [4][4]);
+void get_rtf(struct All_variables *, int, int, double [4][9], int);
+void construct_c3x3matrix_el(struct All_variables *, int, struct CC *, struct CCX *, int, int);
+void get_global_1d_shape_fn(struct All_variables *, int, struct Shape_function1 *, struct Shape_function1_dA *, int);
+void get_global_1d_shape_fn1(struct All_variables *, int, struct Shape_function1 *, struct Shape_function1_dA *, int);
+void mass_matrix(struct All_variables *);
/* Solver_conj_grad.c */
-void set_cg_defaults(struct All_variables *E);
-void cg_allocate_vars(struct All_variables *E);
-void assemble_forces_iterative(struct All_variables *E);
+void cg_allocate_vars(struct All_variables *);
+void assemble_forces_iterative(struct All_variables *);
/* Solver_multigrid.c */
-void set_mg_defaults(struct All_variables *E);
-void mg_allocate_vars(struct All_variables *E);
-void project_vector(struct All_variables *E, int start_lev, double *AU, double *AD, int ic);
-void interp_vector(struct All_variables *E, int start_lev, double *AD, double *AU);
-void project_scalar_e(struct All_variables *E, int start_lev, float *AU, float *AD);
-void project_scalar(struct All_variables *E, int start_lev, float *AU, float *AD);
-void project_viscosity(struct All_variables *E);
-void inject_node_fvector(struct All_variables *E, int start_lev, float **AU, float **AD);
-void inject(struct All_variables *E, int start_lev, double *AU, double *AD);
-void un_inject_vector(struct All_variables *E, int start_lev, double *AD, double *AU);
-void inject_scalar(struct All_variables *E, int start_lev, float *AU, float *AD);
-void inject_scalar_e(struct All_variables *E, int start_lev, float *AU, float *AD);
+void set_mg_defaults(struct All_variables *);
+void mg_allocate_vars(struct All_variables *);
+void project_vector(struct All_variables *, int, double *, double *, int);
+void interp_vector(struct All_variables *, int, double *, double *);
+void project_scalar_e(struct All_variables *, int, float *, float *);
+void project_scalar(struct All_variables *, int, float *, float *);
+void project_viscosity(struct All_variables *);
+void inject_node_fvector(struct All_variables *, int, float **, float **);
+void inject(struct All_variables *, int, double *, double *);
+void un_inject_vector(struct All_variables *, int, double *, double *);
+void inject_scalar(struct All_variables *, int, float *, float *);
+void inject_scalar_e(struct All_variables *, int, float *, float *);
/* Sphere_harmonics.c */
-void set_sphere_harmonics(struct All_variables *E);
-void sphere_harmonics_layer(struct All_variables *E, float **T, float *sphc, float *sphs, int iprint, char *filen);
-void sphere_interpolate(struct All_variables *E, float **T, float *TG);
-void sphere_expansion(struct All_variables *E, float *TG, float *sphc, float *sphs);
-void inv_sphere_harmonics(struct All_variables *E, float *sphc, float *sphs, float *TG, int proc_loc);
-void compute_sphereh_table(struct All_variables *E);
+void sphere_harmonics_layer(struct All_variables *, float **, float *, float *, int, char *);
+void sphere_interpolate(struct All_variables *, float **, float *);
+void sphere_expansion(struct All_variables *, float *, float *, float *);
+void inv_sphere_harmonics(struct All_variables *, float *, float *, float *, int);
+void compute_sphereh_table(struct All_variables *);
/* Stokes_flow_Incomp.c */
-void solve_constrained_flow_iterative(struct All_variables *E);
-float solve_Ahat_p_fhat(struct All_variables *E, double *V, double *P, double *F, double imp, int *steps_max);
-void v_from_vector(struct All_variables *E, float **V, double *F);
+void solve_constrained_flow_iterative(struct All_variables *);
+float solve_Ahat_p_fhat(struct All_variables *, double *, double *, double *, double, int *);
+void v_from_vector(struct All_variables *, float **, double *);
/* Topo_gravity.c */
-void get_CBF_topo(struct All_variables *E, float *H, float *HB);
-void get_STD_topo(struct All_variables *E, float *tpg, float *tpgb, int ii);
+void get_CBF_topo(struct All_variables *, float *, float *);
+void get_STD_topo(struct All_variables *, float *, float *, int);
/* Viscosity_structures.c */
-void viscosity_parameters(struct All_variables *E);
-void get_viscosity_option(struct All_variables *E);
-void viscosity_for_system(struct All_variables *E);
-void get_system_viscosity(struct All_variables *E, int propogate, float *evisc, float *visc);
-void apply_viscosity_smoother(struct All_variables *E, float *visc, float *evisc);
-void visc_from_mat(struct All_variables *E, float *Eta, float *EEta);
-void visc_from_T(struct All_variables *E, float *Eta, float *EEta, int propogate);
-void visc_from_S(struct All_variables *E, float *Eta, float *EEta, int propogate);
-void strain_rate_2_inv(struct All_variables *E, float *EEDOT, int SQRT);
-int layers(struct All_variables *E, float x3);
-int weak_zones(struct All_variables *E, int node, float t_b);
-float boundary_thickness(struct All_variables *E, float *H);
+void viscosity_parameters(struct All_variables *);
+void get_viscosity_option(struct All_variables *);
+void viscosity_for_system(struct All_variables *);
+void get_system_viscosity(struct All_variables *, int, float *, float *);
+void apply_viscosity_smoother(struct All_variables *, float *, float *);
+void visc_from_mat(struct All_variables *, float *, float *);
+void visc_from_T(struct All_variables *, float *, float *, int);
+void visc_from_S(struct All_variables *, float *, float *, int);
+void strain_rate_2_inv(struct All_variables *, float *, int);
+int layers(struct All_variables *, float);
+int weak_zones(struct All_variables *, int, float);
+float boundary_thickness(struct All_variables *, float *);
+void twiddle_thumbs(struct All_variables *, int );
+
+void xyz2rtp(float ,float ,float ,float *);
+void generate_log_message(int ,double ,double ,double , double , struct All_variables *);
+
+void get_ba(struct Shape_function *,struct Shape_function_dx *,
+ struct CC *, struct CCX *, double [4][9],
+ int ,int, double [9][4][9][7]);
+
Modified: mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h 2010-09-15 20:54:02 UTC (rev 17196)
+++ mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h 2010-09-16 05:57:34 UTC (rev 17197)
@@ -51,7 +51,15 @@
int SMOOTH;
int smooth_cycles;
+ int allow_anisotropic_viscosity,anisotropic_viscosity_init;
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ 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
+
char STRUCTURE[20]; /* which option to determine viscosity field, one of .... */
int FROM_SYSTEM;
int FROM_FILE;
@@ -68,6 +76,8 @@
float z410;
float zlith;
float zcomp;
+
+ float zbase_layer[40]; /* new */
int FREEZE;
float freeze_thresh;
More information about the CIG-COMMITS
mailing list