[cig-commits] r17172 - in mc/3D/CitcomS/trunk: . lib
becker at geodynamics.org
becker at geodynamics.org
Sun Sep 5 22:20:44 PDT 2010
Author: becker
Date: 2010-09-05 22:20:44 -0700 (Sun, 05 Sep 2010)
New Revision: 17172
Added:
mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h
Modified:
mc/3D/CitcomS/trunk/configure.ac
mc/3D/CitcomS/trunk/lib/Drive_solvers.c
mc/3D/CitcomS/trunk/lib/Element_calculations.c
mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Makefile.am
mc/3D/CitcomS/trunk/lib/Nodal_mesh.c
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Output_gzdir.c
mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/trunk/lib/Solver_multigrid.c
mc/3D/CitcomS/trunk/lib/Topo_gravity.c
mc/3D/CitcomS/trunk/lib/Tracer_setup.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/ggrd_handling.h
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
Log:
Trial implementation of anisotropic viscosity.
Modified: mc/3D/CitcomS/trunk/configure.ac
===================================================================
--- mc/3D/CitcomS/trunk/configure.ac 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/configure.ac 2010-09-06 05:20:44 UTC (rev 17172)
@@ -103,8 +103,8 @@
# Checks for libraries.
AC_SEARCH_LIBS([MPI_Init], [mpi mpich], [], [AC_MSG_ERROR([MPI library not found])])
-AC_SEARCH_LIBS([sqrt], [m])
+
CIT_CHECK_LIB_HDF5
CIT_CHECK_LIB_HDF5_PARALLEL
@@ -213,6 +213,7 @@
AC_SEARCH_LIBS([ggrd_init_master], [ggrd], [], [AC_MSG_ERROR([HC ggrd library not found. Try setting environment variable HC_HOME.])])
fi
+AC_SEARCH_LIBS([sqrt], [m])
AC_CONFIG_FILES([Makefile
bin/Makefile
Added: mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c (rev 0)
+++ mc/3D/CitcomS/trunk/lib/Anisotropic_viscosity.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -0,0 +1,359 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *<LicenseText>
+ *
+ * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
+ * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
+ * Copyright (C) 1994-2005, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+
+
+/*
+
+ anisotropic viscosity following Muehlhaus, Moresi, Hobbs and Dufour
+ (PAGEOPH, 159, 2311, 2002)
+
+
+
+*/
+
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+
+#include <math.h>
+#include "element_definitions.h"
+#include "global_defs.h"
+#include "material_properties.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))
+/*
+
+
+compute a cartesian anisotropic viscosity matrix
+
+
+ output: D[0,...,5][0,...,5] constitutive matrix
+ input: delta_vis difference in viscosity from isotropic viscosity, which is set to unity
+
+ n[0,..,2]: director orientation, 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],deltac[3][3][3][3];
+
+ /*
+ make sure things are normalized (n[3] might come from interpolation)
+ */
+ nlen = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
+ if((nlen < 0.95)||(nlen > 1.05)){
+ fprintf(stderr,"get_constitutive_orthotropic_viscosity: error: nlen: %g\n",nlen);
+ parallel_process_termination();
+ }
+
+ /* zero out most of D matrix */
+ D[0][1] = D[0][2] = D[0][3] = D[0][4] = D[0][5] = 0.;
+ D[1][0] = D[1][2] = D[1][3] = D[1][4] = D[1][5] = 0.;
+ D[2][0] = D[2][1] = D[2][3] = D[2][4] = D[2][5] = 0.;
+ D[3][0] = D[3][1] = D[3][2] = D[3][4] = D[3][5] = 0.;
+ D[4][0] = D[4][1] = D[4][2] = D[4][3] = D[4][5] = 0.;
+ D[5][0] = D[5][1] = D[5][2] = D[5][3] = D[5][4] = 0.;
+
+ /* isotropic part, in units of iso_visc */
+ 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 */
+
+
+ if(fabs(delta_vis) > 5e-15){
+ /* get Cartesian anisotropy matrix */
+ if(convert_to_spherical){
+ get_delta(deltac,n); /* get anisotropy tensor, \Delta of
+ Muehlhaus et al. (2002) */
+ conv_cart4x4_to_spherical(deltac,theta,phi,delta); /* rotate
+ into
+ CitcomS
+ spherical
+ system */
+ }else{
+ get_delta(delta,n);
+ }
+ delta_vis2 = 2.0*delta_vis;
+ /* s_xx = s_tt */
+ D[0][0] -= delta_vis2 * delta[0][0][0][0]; /* * e_xx */
+ 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]);
+ }
+
+
+}
+
+void set_anisotropic_viscosity_at_element_level(struct All_variables *E, int init_stage)
+{
+ int i,j,k,l,off,nel;
+ double vis2,n[3],u,v,s,r;
+ const int vpts = vpoints[E->mesh.nsd];
+
+ if(E->viscosity.allow_orthotropic_viscosity){
+ if(init_stage){
+ if(E->viscosity.orthotropic_viscosity_init)
+ myerror(E,"anisotropic viscosity should not be initialized twice?!");
+ /* 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.gridmin;i <= E->mesh.gridmax;i++){
+ nel = E->lmesh.NEL[i];
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ 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][j][off] = 0.0;
+ E->EVIn1[i][j][off] = 1.0; E->EVIn2[i][j][off] = E->EVIn3[i][j][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.gridmin;i <= E->mesh.gridmax;i++){
+ nel = E->lmesh.NEL[i];
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ 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][j][off] = vis2;
+ E->EVIn1[i][j][off] = n[0];
+ E->EVIn2[i][j][off] = n[1];
+ E->EVIn3[i][j][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,1);
+ 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.orthotropic_viscosity_init = TRUE;
+ /* end initialization stage */
+ }else{
+ /* standard operation every time step */
+
+
+ }
+ } /* end anisotropic viscosity branch */
+}
+#endif
+
+void normalize_director_at_nodes(struct All_variables *E,float **n1,float **n2, float **n3, int lev)
+{
+ int n,m;
+ double nlen;
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(n=1;n<=E->lmesh.NNO[lev];n++){
+ nlen = sqrt(n1[m][n]*n1[m][n] + n2[m][n]*n2[m][n] + n3[m][n]*n3[m][n]);
+ n1[m][n] /= nlen;
+ n2[m][n] /= nlen;
+ n3[m][n] /= nlen;
+ }
+}
+void normalize_director_at_gint(struct All_variables *E,float **n1,float **n2, float **n3, int lev)
+{
+ int m,e,i,off;
+ const int nsd=E->mesh.nsd;
+ const int vpts=vpoints[nsd];
+ double nlen;
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(e=1;e<=E->lmesh.NEL[lev];e++)
+ for(i=1;i<=vpts;i++) {
+ off = (e-1)*vpts+i;
+ nlen = sqrt(n1[m][off]*n1[m][off] + n2[m][off]*n2[m][off] + n3[m][off]*n3[m][off]);
+ n1[m][off] /= nlen;
+ n2[m][off] /= nlen;
+ n3[m][off] /= nlen;
+ }
+}
+/*
+
+convert cartesian voigt matrix to spherical, CitcomS format
+
+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],base[9];
+ calc_cbase_at_tp_d(theta,phi, base); /* compute cartesian basis at
+ theta, phi location */
+ rot[0][0] = base[3];rot[0][1] = base[4];rot[0][2] = base[5]; /* theta */
+ rot[1][0] = base[6];rot[1][1] = base[7];rot[1][2] = base[8]; /* phi */
+ rot[2][0] = base[0];rot[2][1] = base[1];rot[2][2] = base[2]; /* r */
+ //fprintf(stderr,"%g %g ; %g %g %g ; %g %g %g ; %g %g %g\n\n",theta,phi,rot[0][0],rot[0][1],rot[0][2],rot[1][0],rot[1][1],rot[1][2],rot[2][0],rot[2][1],rot[2][2]);
+ rot_4x4(c,rot,p);
+ //if(E->parallel.me==0)print_6x6_mat(stderr,p);
+ //if(E->parallel.me==0)fprintf(stderr,"\n\n");
+}
+
+/*
+
+
+get anisotropy matrix following Muehlhaus
+
+ */
+void get_delta(double d[3][3][3][3],double n[3])
+{
+ int i,j,k,l;
+ double tmp;
+ 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 */
+void rot_4x4(double c4[3][3][3][3], double r[3][3], double c4c[3][3][3][3])
+{
+
+ int i1,i2,i3,i4,j1,j2,j3,j4;
+ for(i1=0;i1<3;i1++)
+ for(i2=0;i2<3;i2++)
+ for(i3=0;i3<3;i3++)
+ for(i4=0;i4<3;i4++)
+ c4c[i1][i2][i3][i4] = 0.0;
+
+ 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 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");
+ }
+}
+
Modified: mc/3D/CitcomS/trunk/lib/Drive_solvers.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -51,7 +51,6 @@
for (m=1;m<=E->sphere.caps_per_proc;m++)
E->elt_k[i][m]=(struct EK *)malloc((E->lmesh.NEL[i]+1)*sizeof(struct EK));
-
return;
}
@@ -81,7 +80,7 @@
if(need_visc_update(E)){
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);
Modified: mc/3D/CitcomS/trunk/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Element_calculations.c 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Element_calculations.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -46,6 +46,9 @@
static void get_elt_tr(struct All_variables *, int , int , double [24], int );
static void get_elt_tr_pseudo_surf(struct All_variables *, int , int , double [24], int );
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#include "anisotropic_viscosity.h"
+#endif
static void add_force(struct All_variables *E, int e, double elt_f[24], int m)
{
@@ -282,7 +285,7 @@
int lev, iconv;
{
double bdbmu[4][4];
- int pn,qn,ad,bd;
+ int pn,qn,ad,bd,off;
int a,b,i,j,i1,j1,k;
double rtf[4][9],W[9];
@@ -299,76 +302,115 @@
const int ends = ENODES3D;
const int dims=E->mesh.nsd;
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ double D[VPOINTS3D+1][6][6],n[3],btmp[6];
+ int l1,l2;
+#endif
+
get_rtf_at_vpts(E, m, lev, el, rtf);
if (iconv || (el-1)%E->lmesh.ELZ[lev]==0)
construct_c3x3matrix_el(E,el,&E->element_Cc,&E->element_Ccx,lev,m,0);
-
+
/* Note N[a].gauss_pt[n] is the value of shape fn a at the nth gaussian
quadrature point. Nx[d] is the derivative wrt x[d]. */
-
+
for(k=1;k<=vpts;k++) {
- W[k]=g_point[k].weight[dims-1]*E->GDA[lev][m][el].vpt[k]*E->EVI[lev][m][(el-1)*vpts+k];
+ off = (el-1)*vpts+k;
+ W[k]=g_point[k].weight[dims-1]*E->GDA[lev][m][el].vpt[k]*E->EVI[lev][m][off];
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ if(E->viscosity.allow_orthotropic_viscosity){/* allow for a possibly anisotropic viscosity */
+ n[0] = E->EVIn1[lev][m][off];n[1] = E->EVIn2[lev][m][off];n[2] = E->EVIn3[lev][m][off]; /* Cartesian directors */
+ /* get the viscosity factor matrix and convert to CitcomS spherical */
+ get_constitutive_orthotropic_viscosity(D[k],E->EVI2[lev][m][off],n,TRUE,rtf[1][k],rtf[2][k]);
+ }
+#endif
}
-
+
get_ba(&(E->N), &(E->GNX[lev][m][el]), &E->element_Cc, &E->element_Ccx,
rtf, E->mesh.nsd, ba);
- 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;
+ for(a=1;a<=ends;a++) /* loop over element nodes */
+ for(b=a;b<=ends;b++) {
- 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][k][i][1]*ba[b][k][j][1] +
- ba[a][k][i][2]*ba[b][k][j][2] +
- ba[a][k][i][3]*ba[b][k][j][3] ) +
- ba[a][k][i][4]*ba[b][k][j][4] +
- ba[a][k][i][5]*ba[b][k][j][5] +
- ba[a][k][i][6]*ba[b][k][j][6] );
+ 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;
- if(E->control.inv_gruneisen != 0)
- for(i=1;i<=dims;i++)
- for(j=1;j<=dims;j++)
- for(k=1;k<=VPOINTS3D;k++)
- bdbmu[i][j] -= W[k] * two_thirds *
- ( ba[a][k][i][1] + ba[a][k][i][2] + ba[a][k][i][3] ) *
- ( ba[b][k][j][1] + ba[b][k][j][2] + ba[b][k][j][3] );
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ if(E->viscosity.allow_orthotropic_viscosity){
+ for(i=1;i<=dims;i++)
+ for(j=1;j<=dims;j++)
+ for(k=1;k<=vpts;k++){ /* */
-
- /**/
- ad=dims*(a-1);
- bd=dims*(b-1);
-
- pn=ad*nn+bd;
- qn=bd*nn+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+nn ] = bdbmu[2][1] ;
- elt_k[pn+nn+1 ] = bdbmu[2][2] ;
- elt_k[pn+nn+2 ] = bdbmu[2][3] ;
- elt_k[pn+2*nn ] = bdbmu[3][1] ;
- elt_k[pn+2*nn+1] = bdbmu[3][2] ;
- elt_k[pn+2*nn+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+nn ] = bdbmu[1][2] ;
- elt_k[qn+nn+1 ] = bdbmu[2][2] ;
- elt_k[qn+nn+2 ] = bdbmu[3][2] ;
- elt_k[qn+2*nn ] = bdbmu[1][3] ;
- elt_k[qn+2*nn+1] = bdbmu[2][3] ;
- elt_k[qn+2*nn+2] = bdbmu[3][3] ;
- /**/
-
+ /* note that D is in 0,...,N-1 convention */
+ for(l1=0;l1 < 6;l1++) /* compute D*B */
+ for(btmp[l1]=0.0,l2=0;l2<6;l2++)
+ btmp[l1] += D[k][l1][l2] * ba[b][k][j][l2+1];
+ /* compute B^T (D*B) */
+ bdbmu[i][j] += W[k] * ( ba[a][k][i][1]*btmp[0]+ba[a][k][i][2]*btmp[1]+ba[a][k][i][3]*btmp[2]+
+ ba[a][k][i][4]*btmp[3]+ba[a][k][i][5]*btmp[4]+ba[a][k][i][6]*btmp[5]);
+ }
+ if(E->control.inv_gruneisen != 0)
+ for(i=1;i<=dims;i++)
+ for(j=1;j<=dims;j++)
+ for(k=1;k<=VPOINTS3D;k++)
+ bdbmu[i][j] -= W[k] * two_thirds *
+ ( ba[a][k][i][1] + ba[a][k][i][2] + ba[a][k][i][3] ) *
+ ( ba[b][k][j][1] + ba[b][k][j][2] + ba[b][k][j][3] );
+ }else{
+#endif /* isotropic branch */
+ 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][k][i][1]*ba[b][k][j][1] +
+ ba[a][k][i][2]*ba[b][k][j][2] +
+ ba[a][k][i][3]*ba[b][k][j][3] ) +
+ ba[a][k][i][4]*ba[b][k][j][4] +
+ ba[a][k][i][5]*ba[b][k][j][5] +
+ ba[a][k][i][6]*ba[b][k][j][6] );
+
+ if(E->control.inv_gruneisen != 0)
+ for(i=1;i<=dims;i++)
+ for(j=1;j<=dims;j++)
+ for(k=1;k<=VPOINTS3D;k++)
+ bdbmu[i][j] -= W[k] * two_thirds *
+ ( ba[a][k][i][1] + ba[a][k][i][2] + ba[a][k][i][3] ) *
+ ( ba[b][k][j][1] + ba[b][k][j][2] + ba[b][k][j][3] );
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ }
+#endif
+
+ /**/
+ ad=dims*(a-1);
+ bd=dims*(b-1);
+
+ pn=ad*nn+bd;
+ qn=bd*nn+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+nn ] = bdbmu[2][1] ;
+ elt_k[pn+nn+1 ] = bdbmu[2][2] ;
+ elt_k[pn+nn+2 ] = bdbmu[2][3] ;
+ elt_k[pn+2*nn ] = bdbmu[3][1] ;
+ elt_k[pn+2*nn+1] = bdbmu[3][2] ;
+ elt_k[pn+2*nn+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+nn ] = bdbmu[1][2] ;
+ elt_k[qn+nn+1 ] = bdbmu[2][2] ;
+ elt_k[qn+nn+2 ] = bdbmu[3][2] ;
+ elt_k[qn+2*nn ] = bdbmu[1][3] ;
+ elt_k[qn+2*nn+1] = bdbmu[2][3] ;
+ elt_k[qn+2*nn+2] = bdbmu[3][3] ;
+ /**/
+
} /* Sum over all the a,b's to obtain full elt_k matrix */
-
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -44,6 +44,9 @@
#include "composition_related.h"
#include "element_definitions.h"
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#include "anisotropic_viscosity.h"
+#endif
#ifdef USE_GGRD
#include "hc.h" /* ggrd and hc packages */
@@ -56,6 +59,9 @@
int layers(struct All_variables *,int ,int );
void ggrd_vtop_helper_decide_on_internal_nodes(struct All_variables *, /* input */
int ,int ,int, int,int ,int *, int *,int *);
+void convert_pvec_to_cvec_d(double ,double , double , double *,double *);
+void calc_cbase_at_tp_d(double , double , double *);
+void xyz2rtpd(float ,float ,float ,double *);
/*
@@ -132,11 +138,13 @@
phi*180/M_PI, 90-theta*180/M_PI);
myerror(E,error);
}
- /* limit to 0 or 1 */
- if(indbl < .5)
- indbl = 0.0;
- else
- indbl = 1.0;
+ if(!E->control.ggrd_comp_smooth){
+ /* limit to 0 or 1 */
+ if(indbl < .5)
+ indbl = 0.0;
+ else
+ indbl = 1.0;
+ }
E->trace.extraq[j][0][kk]= indbl;
}else{
/* below */
@@ -359,7 +367,7 @@
int mpi_rc,timedep,interpolate;
int mpi_inmsg, mpi_success_message = 1;
int m,el,i,j,k,inode,i1,i2,elxlz,elxlylz,ind;
- int llayer,nox,noy,noz,nox1,noz1,noy1,level,lselect,idim,elx,ely,elz;
+ int llayer,nox,noy,noz,level,lselect,idim,elx,ely,elz;
char gmt_string[10],char_dummy;
double indbl,indbl2,age,f1,f2,vip,rout[3],xloc[4];
char tfilename[1000];
@@ -369,7 +377,6 @@
FILE *in;
nox=E->mesh.nox;noy=E->mesh.noy;noz=E->mesh.noz;
- nox1=E->lmesh.nox;noz1=E->lmesh.noz;noy1=E->lmesh.noy;
elx=E->lmesh.elx;elz=E->lmesh.elz;ely=E->lmesh.ely;
elxlz = elx * elz;
elxlylz = elxlz * ely;
@@ -577,7 +584,7 @@
int mpi_rc,timedep,interpolate;
int mpi_inmsg, mpi_success_message = 1;
int m,el,i,j,k,node,i1,i2,elxlz,elxlylz,ind;
- int llayer,nox,noy,noz,nox1,noz1,noy1,lev,lselect,idim,elx,ely,elz;
+ int llayer,nox,noy,noz,lev,lselect,idim,elx,ely,elz;
char gmt_string[10],char_dummy;
double indbl,indbl2,age,f1,f2,vip,rout[3],xloc[4];
char tfilename[1000];
@@ -587,7 +594,6 @@
const int ends = enodes[dims];
/* dimensional ints */
nox=E->mesh.nox;noy=E->mesh.noy;noz=E->mesh.noz;
- nox1=E->lmesh.nox;noz1=E->lmesh.noz;noy1=E->lmesh.noy;
elx=E->lmesh.elx;elz=E->lmesh.elz;ely=E->lmesh.ely;
elxlz = elx * elz;
elxlylz = elxlz * ely;
@@ -1306,6 +1312,210 @@
+/*
+
+
+read in anisotropic viscosity from a directory which holds
+
+
+vis2.grd for the viscosity factors
+
+nr.grd, nt.grd, np.grd for the directors
+
+*/
+void ggrd_read_anivisc_from_file(struct All_variables *E, int is_global)
+{
+ MPI_Status mpi_stat;
+ int mpi_rc;
+ int mpi_inmsg, mpi_success_message = 1;
+ int m,el,i,j,k,l,inode,i1,i2,elxlz,elxlylz,ind,nel;
+ int llayer,nox,noy,noz,level,lselect,idim,elx,ely,elz;
+ char gmt_string[10],char_dummy;
+ double vis2,ntheta,nphi,nr,rout[3],xloc[4],nlen;
+ double cvec[3],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 init_layer = 1; /* initialize all elements in top layer */
+ 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;
+
+#ifndef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ fprintf(stderr,"ggrd_read_anivisc_from_file: error, need to compile with CITCOM_ALLOW_ORTHOTROPIC_VISC\n");
+ parallel_process_termination();
#endif
+ if(!E->viscosity.allow_orthotropic_viscosity)
+ myerror(E,"ggrd_read_anivisc_from_file: called, but allow_orthotropic_viscosity is FALSE?!");
+
+ /* isotropic default */
+ for(i=E->mesh.gridmin;i <= E->mesh.gridmax;i++){
+ nel = E->lmesh.NEL[i];
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ 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][j][ind] = 0.0;
+ E->EVIn1[i][j][ind] = 1.0; E->EVIn2[i][j][ind] = E->EVIn3[i][j][ind] = 0.0;
+ }
+ }
+ }
+ }
+ /*
+
+ rest
+ */
+ if(is_global) /* decide on GMT flag */
+ sprintf(gmt_string,GGRD_GMT_GLOBAL_STRING); /* global */
+ else
+ sprintf(gmt_string,"");
+ /*
+
+ 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, E->parallel.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)
+ fprintf(stderr,"ggrd_read_anivisc_from_file: initializing, assigning to all above %g km, input is %s\n",
+ E->data.radius_km*E->viscosity.zbase_layer[E->control.ggrd.mat_control-1],
+ (is_global)?("global"):("regional"));
+ /*
+ 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(E,"ggrd init error");
+ /* n_r */
+ 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(E,"ggrd init error");
+ /* n_theta */
+ 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(E,"ggrd init error");
+ /* n_phi */
+ 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(E,"ggrd init error");
+
+ /* 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, E->parallel.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 (m=1;m <= E->sphere.caps_per_proc;m++) {
+ for (j=1;j <= elz;j++) { /* this assumes a regular grid sorted as in (1)!!! */
+ if(E->mat[m][j] <= init_layer){
+ /* 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
+ */
+ xloc[1] = xloc[2] = xloc[3] = 0.0;
+ for(inode=1;inode <= 4;inode++){
+ ind = E->ien[m][el].node[inode];
+ xloc[1] += E->x[m][1][ind];xloc[2] += E->x[m][2][ind];xloc[3] += E->x[m][3][ind];
+ }
+ xloc[1]/=4.;xloc[2]/=4.;xloc[3]/=4.;
+ xyz2rtpd(xloc[1],xloc[2],xloc[3],rout); /* convert to spherical */
+
+ /* vis2 */
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],
+ vis2_grd,&vis2,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();
+ }
+ /*
+ convert from n_east,n_north to Cartesian vector,
+ assuming the director is in the horizontal,
+ i.e. n_r = 0
+ */
+ nlen = sqrt(nphi*nphi+ntheta*ntheta+nr*nr); /* correct, because
+ interpolation might have
+ screwed up initialization */
+ nphi /= nlen; ntheta /= nlen;nr /= nlen;
+ calc_cbase_at_tp_d(rout[1],rout[2],base);
+ convert_pvec_to_cvec_d(nr,ntheta,nphi,base,cvec);
+ for(l=1;l <= vpts;l++){ /* assign to all integration points */
+ ind = (el-1)*vpts + l;
+ E->EVI2[E->mesh.gridmax][m][ind] = vis2;
+ E->EVIn1[E->mesh.gridmax][m][ind] = cvec[0];
+ E->EVIn2[E->mesh.gridmax][m][ind] = cvec[1];
+ E->EVIn3[E->mesh.gridmax][m][ind] = cvec[2];
+ }
+ }
+ }
+ } /* end insize lith */
+ } /* end elz loop */
+ } /* end m 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
+
+
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -160,7 +160,11 @@
construct_sub_element(E);
construct_shape_functions(E);
construct_shape_function_derivatives(E);
+ if(chatty)fprintf(stderr,"shape functions done\n");
+
construct_elt_gs(E);
+
+
if(E->control.inv_gruneisen != 0)
construct_elt_cs(E);
@@ -589,7 +593,12 @@
*/
input_boolean("allow_mixed_vbcs",&(E->control.ggrd_allow_mixed_vbcs),"off",m);
+ /* when assigning composition from grid file, allow values between 0 and 1 ?
+ default will ronud up/down to 0 or 1
+ */
+ input_boolean("ggrd_comp_smooth",&(E->control.ggrd_comp_smooth),"off",m);
+
#endif
input_boolean("aug_lagr",&(E->control.augmented_Lagr),"off",m);
@@ -878,7 +887,7 @@
{
void set_up_nonmg_aliases();
- int m,n,snel,nsf,elx,ely,nox,noy,noz,nno,nel,npno;
+ int m,n,snel,nsf,elx,ely,nox,noy,noz,nno,nel,npno,lim;
int k,i,j,d,l,nno_l,npno_l,nozl,nnov_l,nxyz;
m=0;
@@ -1045,6 +1054,34 @@
} /* end for cap and i & j */
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ if(E->viscosity.allow_orthotropic_viscosity){
+ for(i=E->mesh.gridmin;i<=E->mesh.gridmax;i++)
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ nel = E->lmesh.NEL[i];
+ nno = E->lmesh.NNO[i];
+ E->EVI2[i][j] = (float *) malloc((nel+1)*vpoints[E->mesh.nsd]*sizeof(float));
+ E->EVIn1[i][j] = (float *) malloc((nel+1)*vpoints[E->mesh.nsd]*sizeof(float));
+ E->EVIn2[i][j] = (float *) malloc((nel+1)*vpoints[E->mesh.nsd]*sizeof(float));
+ E->EVIn3[i][j] = (float *) malloc((nel+1)*vpoints[E->mesh.nsd]*sizeof(float));
+
+ E->VI2[i][j] = (float *) malloc((nno+1)*sizeof(float));
+ E->VIn1[i][j] = (float *) malloc((nno+1)*sizeof(float));
+ E->VIn2[i][j] = (float *) malloc((nno+1)*sizeof(float));
+ E->VIn3[i][j] = (float *) malloc((nno+1)*sizeof(float));
+ if((!(E->EVI2[i][j]))||(!(E->VI2[i][j]))||
+ (!(E->EVIn1[i][j]))||(!(E->EVIn2[i][j]))||(!(E->EVIn3[i][j]))||
+ (!(E->VIn1[i][j]))||(!(E->VIn2[i][j]))||(!(E->VIn3[i][j]))){
+ fprintf(stderr, "Error: Cannot allocate anisotropic visc memory, rank=%d\n",
+ E->parallel.me);
+ parallel_process_termination();
+ }
+ }
+ E->viscosity.orthotropic_viscosity_init = FALSE;
+ if(E->parallel.me == 0)
+ fprintf(stderr,"allocated for anisotropic viscosity levmax %i\n",E->mesh.gridmax);
+ }
+#endif
for (j=1;j<=E->sphere.caps_per_proc;j++) {
Modified: mc/3D/CitcomS/trunk/lib/Makefile.am
===================================================================
--- mc/3D/CitcomS/trunk/lib/Makefile.am 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Makefile.am 2010-09-06 05:20:44 UTC (rev 17172)
@@ -68,7 +68,9 @@
sources = \
Advection_diffusion.c \
+ Anisotropic_viscosity.c \
advection_diffusion.h \
+ anisotropic_viscosity.h \
advection.h \
BC_util.c \
Checkpoints.c \
Modified: mc/3D/CitcomS/trunk/lib/Nodal_mesh.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Nodal_mesh.c 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Nodal_mesh.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -34,6 +34,7 @@
#include "global_defs.h"
+
/* get nodal spherical velocities from the solution vector */
void v_from_vector(E)
struct All_variables *E;
@@ -230,127 +231,132 @@
}
+/*
+ interpolate the viscosity from element integration points to nodes
+
+ */
void visc_from_gint_to_nodes(E,VE,VN,lev)
struct All_variables *E;
float **VE,**VN;
int lev;
- {
- int m,e,i,j,k,n;
+{
+ int m,e,i,j,k,n,off,lim;
const int nsd=E->mesh.nsd;
const int vpts=vpoints[nsd];
const int ends=enodes[nsd];
- double temp_visc;
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.NNO[lev];i++)
- VN[m][i] = 0.0;
-
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(e=1;e<=E->lmesh.NEL[lev];e++) {
- temp_visc=0.0;
- for(i=1;i<=vpts;i++)
- temp_visc += VE[m][(e-1)*vpts + i];
- temp_visc = temp_visc/vpts;
-
- for(j=1;j<=ends;j++) {
- n = E->IEN[lev][m][e].node[j];
- VN[m][n] += E->TWW[lev][m][e].node[j] * temp_visc;
- }
+ double temp_visc, weight;
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=E->lmesh.NNO[lev];i++)
+ VN[m][i] = 0.0;
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(e=1;e<=E->lmesh.NEL[lev];e++) {
+ temp_visc=0.0;
+ for(i=1;i<=vpts;i++)
+ temp_visc += VE[m][(e-1)*vpts + i];
+ temp_visc = temp_visc/vpts;
+
+ for(j=1;j<=ends;j++) {
+ n = E->IEN[lev][m][e].node[j];
+ VN[m][n] += E->TWW[lev][m][e].node[j] * temp_visc;
+ }
}
+ (E->exchange_node_f)(E,VN,lev);
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(n=1;n<=E->lmesh.NNO[lev];n++)
+ VN[m][n] *= E->MASS[lev][m][n];
- (E->exchange_node_f)(E,VN,lev);
+ return;
+}
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(n=1;n<=E->lmesh.NNO[lev];n++)
- VN[m][n] *= E->MASS[lev][m][n];
+/*
- return;
-}
+interpolate viscosity from nodes to element integration points
-
+ */
void visc_from_nodes_to_gint(E,VN,VE,lev)
struct All_variables *E;
float **VE,**VN;
int lev;
- {
+{
- int m,e,i,j,k,n;
+ int m,e,i,j,k,n,off;
const int nsd=E->mesh.nsd;
const int vpts=vpoints[nsd];
const int ends=enodes[nsd];
double temp_visc;
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(e=1;e<=E->lmesh.NEL[lev];e++)
- for(i=1;i<=vpts;i++)
- VE[m][(e-1)*vpts+i] = 0.0;
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(e=1;e<=E->lmesh.NEL[lev];e++)
- for(i=1;i<=vpts;i++) {
- temp_visc=0.0;
- for(j=1;j<=ends;j++)
- temp_visc += E->N.vpt[GNVINDEX(j,i)]*VN[m][E->IEN[lev][m][e].node[j]];
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(e=1;e<=E->lmesh.NEL[lev];e++)
+ for(i=1;i<=vpts;i++)
+ VE[m][(e-1)*vpts+i] = 0.0;
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(e=1;e<=E->lmesh.NEL[lev];e++)
+ for(i=1;i<=vpts;i++) {
+ temp_visc=0.0;
+ for(j=1;j<=ends;j++)
+ temp_visc += E->N.vpt[GNVINDEX(j,i)]*VN[m][E->IEN[lev][m][e].node[j]];
+ VE[m][(e-1)*vpts+i] = temp_visc;
+ }
- VE[m][(e-1)*vpts+i] = temp_visc;
- }
+ return;
+}
- return;
- }
+/* called from MG as (?)
+ visc_from_gint_to_ele(E,E->EVI[lv],viscU,lv)
+
+*/
void visc_from_gint_to_ele(E,VE,VN,lev)
struct All_variables *E;
float **VE,**VN;
int lev;
{
- int m,e,i,j,k,n;
- const int nsd=E->mesh.nsd;
- const int vpts=vpoints[nsd];
- const int ends=enodes[nsd];
- double temp_visc;
+ int m,e,i,j,k,n,off;
+ const int nsd=E->mesh.nsd;
+ const int vpts=vpoints[nsd];
+ const int ends=enodes[nsd];
+ double temp_visc;
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.NEL[lev];i++)
- VN[m][i] = 0.0;
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=E->lmesh.NEL[lev];i++)
+ VN[m][i] = 0.0;
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(e=1;e<=E->lmesh.NEL[lev];e++) {
+ temp_visc=0.0;
+ for(i=1;i<=vpts;i++)
+ temp_visc += VE[m][(e-1)*vpts + i];
+ temp_visc = temp_visc/vpts;
+ VN[m][e] = temp_visc;
+ }
+
+ return;
+ }
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(e=1;e<=E->lmesh.NEL[lev];e++) {
- temp_visc=0.0;
- for(i=1;i<=vpts;i++)
- temp_visc += VE[m][(e-1)*vpts + i];
- temp_visc = temp_visc/vpts;
+/* called from MG as
- VN[m][e] = temp_visc;
- }
+ visc_from_ele_to_gint(E,viscD,E->EVI[sl_minus],sl_minus);
- return;
-}
+*/
-
void visc_from_ele_to_gint(E,VN,VE,lev)
struct All_variables *E;
float **VE,**VN;
int lev;
- {
-
- int m,e,i,j,k,n;
+{
+ int m,e,i,j,k,n,off;
const int nsd=E->mesh.nsd;
const int vpts=vpoints[nsd];
const int ends=enodes[nsd];
double temp_visc;
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(e=1;e<=E->lmesh.NEL[lev];e++)
- for(i=1;i<=vpts;i++)
- VE[m][(e-1)*vpts+i] = 0.0;
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(e=1;e<=E->lmesh.NEL[lev];e++)
- for(i=1;i<=vpts;i++) {
-
- VE[m][(e-1)*vpts+i] = VN[m][e];
- }
-
- return;
- }
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(e=1;e<=E->lmesh.NEL[lev];e++)
+ for(i=1;i<=vpts;i++) {
+ VE[m][(e-1)*vpts+i] = VN[m][e];
+ }
+ return;
+}
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -58,7 +58,14 @@
extern void get_STD_topo(struct All_variables *, float**, float**,
float**, float**, int);
extern void get_CBF_topo(struct All_variables *, float**, float**);
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#include "anisotropic_viscosity.h"
+void output_avisc(struct All_variables *, int);
+#endif
+
+
+
/**********************************************************************/
void output_common_input(struct All_variables *E)
@@ -105,7 +112,11 @@
output_velo(E, cycles);
output_visc(E, cycles);
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ output_avisc(E, cycles);
+#endif
+
output_surf_botm(E, cycles);
/* optional output below */
@@ -311,7 +322,30 @@
return;
}
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+void output_avisc(struct All_variables *E, int cycles)
+{
+ int i, j;
+ char output_file[255];
+ FILE *fp1;
+ int lev = E->mesh.levmax;
+ if(E->viscosity.allow_orthotropic_viscosity){
+ sprintf(output_file,"%s.avisc.%d.%d", E->control.data_file,
+ E->parallel.me, cycles);
+ fp1 = output_open(output_file, "w");
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ fprintf(fp1,"%3d %7d\n",j,E->lmesh.nno);
+ for(i=1;i<=E->lmesh.nno;i++)
+ fprintf(fp1,"%.4e %.4e %.4e %.4e\n",E->VI2[lev][j][i],E->VIn1[lev][j][i],E->VIn2[lev][j][i],E->VIn3[lev][j][i]);
+ }
+
+ fclose(fp1);
+ }
+ return;
+}
+#endif
+
void output_velo(struct All_variables *E, int cycles)
{
int i, j;
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -117,6 +117,10 @@
int open_file_zipped(char *, FILE **,struct All_variables *);
void gzip_file(char *);
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#include "anisotropic_viscosity.h"
+void gzdir_output_avisc(struct All_variables *, int);
+#endif
extern void temperatures_conform_bcs(struct All_variables *);
extern void myerror(struct All_variables *,char *);
@@ -163,6 +167,9 @@
else new VTK output won't
work */
gzdir_output_visc(E, out_cycles);
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ gzdir_output_avisc(E, out_cycles);
+#endif
gzdir_output_surf_botm(E, out_cycles);
@@ -753,8 +760,68 @@
return;
}
+/*
+ anisotropic viscosity
+*/
+void gzdir_output_avisc(struct All_variables *E, int cycles)
+{
+ int i, j;
+ char output_file[255];
+ gzFile *gz1;
+ FILE *fp1;
+ int lev = E->mesh.levmax;
+ float ftmp;
+ /* for dealing with several processors */
+ MPI_Status mpi_stat;
+ int mpi_rc;
+ int mpi_inmsg, mpi_success_message = 1;
+ if(E->viscosity.allow_orthotropic_viscosity){
+
+ if(E->output.gzdir.vtk_io < 2){
+ snprintf(output_file,255,
+ "%s/%d/avisc.%d.%d.gz", E->control.data_dir,
+ cycles,E->parallel.me, cycles);
+ gz1 = gzdir_output_open(output_file,"w");
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ gzprintf(gz1,"%3d %7d\n",j,E->lmesh.nno);
+ for(i=1;i<=E->lmesh.nno;i++)
+ gzprintf(gz1,"%.4e %.4e %.4e %.4e\n",E->VI2[lev][j][i],E->VIn1[lev][j][i],E->VIn2[lev][j][i],E->VIn3[lev][j][i]);
+ }
+
+ gzclose(gz1);
+ }else{
+ if(E->output.gzdir.vtk_io == 2)
+ parallel_process_sync(E);
+ /* new legacy VTK */
+ get_vtk_filename(output_file,0,E,cycles);
+ if((E->parallel.me == 0) || (E->output.gzdir.vtk_io == 3)){
+ fp1 = output_open(output_file,"a");
+ myfprintf(fp1,"SCALARS vis2 float 1\n");
+ myfprintf(fp1,"LOOKUP_TABLE default\n");
+ }else{
+ /* if not first, wait for previous */
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), 0, E->parallel.world, &mpi_stat);
+ /* open for append */
+ fp1 = output_open(output_file,"a");
+ }
+ for(j=1; j<= E->sphere.caps_per_proc;j++)
+ for(i=1;i<=E->lmesh.nno;i++){
+ ftmp = E->VI2[lev][j][i];
+ if(be_write_float_to_file(&ftmp,1,fp1)!=1)
+ BE_WERROR;
+ }
+ fclose(fp1);fflush(fp1); /* close file and flush buffer */
+ if(E->output.gzdir.vtk_io == 2)
+ if(E->parallel.me < E->parallel.nproc-1){
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 0, E->parallel.world);
+ }
+ }
+ }
+ return;
+}
+
void gzdir_output_surf_botm(struct All_variables *E, int cycles)
{
int i, j, s;
Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -52,9 +52,11 @@
#include "parallel_related.h"
void calc_cbase_at_tp(float , float , float *);
+void calc_cbase_at_tp_d(double , double , double *);
void rtp2xyz(float , float , float, float *);
void rtp2xyzd(double , double , double, double *);
void convert_pvec_to_cvec(float ,float , float , float *,float *);
+void convert_pvec_to_cvec_d(double ,double , double , double *,double *);
void *safe_malloc (size_t );
void myerror(struct All_variables *,char *);
void xyz2rtp(float ,float ,float ,float *);
@@ -423,7 +425,30 @@
base[7]= cp;
base[8]= 0.0;
}
+void calc_cbase_at_tp_d(double theta, double phi, double *base) /* double version */
+{
+
+ double ct,cp,st,sp;
+
+ ct=cos(theta);
+ cp=cos(phi);
+ st=sin(theta);
+ sp=sin(phi);
+ /* r */
+ base[0]= st * cp;
+ base[1]= st * sp;
+ base[2]= ct;
+ /* theta */
+ base[3]= ct * cp;
+ base[4]= ct * sp;
+ base[5]= -st;
+ /* phi */
+ base[6]= -sp;
+ base[7]= cp;
+ base[8]= 0.0;
+}
+
/* calculate base at nodal locations where we have precomputed cos/sin */
void calc_cbase_at_node(int cap, int node, float *base,struct All_variables *E)
@@ -463,6 +488,17 @@
cvec[i] += base[6+i]* vp;
}
}
+void convert_pvec_to_cvec_d(double vr,double vt,
+ double vp, double *base,
+ double *cvec)
+{
+ int i;
+ for(i=0;i<3;i++){
+ cvec[i] = base[i] * vr;
+ cvec[i] += base[3+i]* vt;
+ cvec[i] += base[6+i]* vp;
+ }
+}
/*
like malloc, but with test
Modified: mc/3D/CitcomS/trunk/lib/Solver_multigrid.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Solver_multigrid.c 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Solver_multigrid.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -253,32 +253,82 @@
viscD[m]=(float *)malloc((1+vpts*E->lmesh.NEL[lv-1])*sizeof(float));
}
- 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);
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC /* allow for anisotropy */
+ if(E->viscosity.allow_orthotropic_viscosity){
+ 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); /* isotropic */
+ project_scalar(E,lv,viscU,viscD);
+ visc_from_nodes_to_gint(E,viscD,E->EVI[sl_minus],sl_minus);
+ /* anisotropic */
+ 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);
+ normalize_director_at_gint(E,E->EVIn1[sl_minus],E->EVIn2[sl_minus],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]);
+ else if (E->viscosity.smooth_cycles==2) {
+ visc_from_gint_to_ele(E,E->EVI[lv],viscU,lv); /* isotropic */
+ inject_scalar_e(E,lv,viscU,E->EVI[sl_minus]);
+ /* anisotropic */
+ 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]);
+ normalize_director_at_gint(E,E->EVIn1[sl_minus],E->EVIn2[sl_minus],E->EVIn3[sl_minus],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);
+ else if (E->viscosity.smooth_cycles==3) {
+ visc_from_gint_to_ele(E,E->EVI[lv],viscU,lv);/* isotropic */
+ project_scalar_e(E,lv,viscU,viscD);
+ visc_from_ele_to_gint(E,viscD,E->EVI[sl_minus],sl_minus);
+ /* anisotropic */
+ 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);
+ normalize_director_at_gint(E,E->EVIn1[sl_minus],E->EVIn2[sl_minus],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); */
-
- inject_scalar(E,lv,E->VI[lv],E->VI[sl_minus]);
- visc_from_nodes_to_gint(E,E->VI[sl_minus],E->EVI[sl_minus],sl_minus);
+ else if (E->viscosity.smooth_cycles==0) {
+ inject_scalar(E,lv,E->VI[lv],E->VI[sl_minus]);/* isotropic */
+ visc_from_nodes_to_gint(E,E->VI[sl_minus],E->EVI[sl_minus],sl_minus);
+ /* anisotropic */
+ inject_scalar(E,lv,E->VI2[lv],E->VI2[sl_minus]);visc_from_nodes_to_gint(E,E->VI2[sl_minus],E->EVI2[sl_minus],sl_minus);
+ inject_scalar(E,lv,E->VIn1[lv],E->VIn1[sl_minus]);visc_from_nodes_to_gint(E,E->VIn1[sl_minus],E->EVIn1[sl_minus],sl_minus);
+ inject_scalar(E,lv,E->VIn2[lv],E->VIn2[sl_minus]);visc_from_nodes_to_gint(E,E->VIn2[sl_minus],E->EVIn2[sl_minus],sl_minus);
+ inject_scalar(E,lv,E->VIn3[lv],E->VIn3[sl_minus]);visc_from_nodes_to_gint(E,E->VIn3[sl_minus],E->EVIn3[sl_minus],sl_minus);
+ normalize_director_at_gint(E,E->EVIn1[sl_minus],E->EVIn2[sl_minus],E->EVIn3[sl_minus],sl_minus);
}
+ }
+ }else{
+#endif
+ /* regular operation, isotropic viscosity */
+ 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);
+ }
+ 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]);
+ }
+ 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);
+ }
+ 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); */
+
+ inject_scalar(E,lv,E->VI[lv],E->VI[sl_minus]);
+ visc_from_nodes_to_gint(E,E->VI[sl_minus],E->EVI[sl_minus],sl_minus);
+ }
/* for(m=1;m<=E->sphere.caps_per_proc;m++) {
@@ -289,7 +339,9 @@
}
*/
}
-
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ }
+#endif
for(m=1;m<=E->sphere.caps_per_proc;m++) {
free((void *)viscU[m]);
free((void *)viscD[m]);
Modified: mc/3D/CitcomS/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -48,7 +48,11 @@
float** , float** , float** ,
float** , float** );
void stress_conform_bcs(struct All_variables *);
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#include "anisotropic_viscosity.h"
+#endif
+
/*
compute the full stress tensor and the dynamic topo
@@ -217,7 +221,7 @@
void construct_c3x3matrix_el();
void get_ba();
- int i,j,p,q,e,node,m;
+ int i,j,p,q,e,node,m,l1,l2;
float VV[4][9],Vxyz[9][9],Szz,Sxx,Syy,Sxy,Sxz,Szy,div,vor;
double dilation[9];
@@ -225,7 +229,9 @@
double pre[9],tww[9],rtf[4][9];
double velo_scaling, stress_scaling, mass_fac;
-
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ double D[6][6],n[3],eps[6],str[6];
+#endif
struct Shape_function_dA *dOmega;
struct Shape_function_dx *GNx;
@@ -269,7 +275,7 @@
* where 1 is theta, 2 is phi, and 3 is r
*/
for(j=1;j <= vpts;j++) { /* loop through velocity Gauss points */
-
+ /* E->EVi[j] = E->EVI[E->mesh.levmax][j]; */
pre[j] = E->EVi[m][(e-1)*vpts+j]*dOmega->vpt[j];
dilation[j] = 0.0;
Vxyz[1][j] = 0.0;
@@ -359,7 +365,43 @@
for(j=1;j<=vpts;j++)
dilation[j] = (Vxyz[1][j] + Vxyz[2][j] + Vxyz[3][j]) / 3.0;
}
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ if(E->viscosity.allow_orthotropic_viscosity){
+ for(j=1;j <= vpts;j++) {
+ l1 = (e-1)*vpts+j;
+ n[0] = E->EVIn1[E->mesh.levmax][m][l1]; /* Cartesian directors */
+ n[1] = E->EVIn2[E->mesh.levmax][m][l1];
+ n[2] = E->EVIn3[E->mesh.levmax][m][l1];
+ /*
+ get viscosity matrix and convert to spherical system in
+ CitcomS convection
+ */
+ get_constitutive_orthotropic_viscosity(D,E->EVI2[E->mesh.levmax][m][l1],n,TRUE,rtf[1][j],rtf[2][j]);
+
+ /* deviatoric stress, pressure will be added later */
+ eps[0] = Vxyz[1][j] - dilation[j]; /* strain-rates */
+ eps[1] = Vxyz[2][j] - dilation[j];
+ eps[2] = Vxyz[3][j] - dilation[j];
+ eps[3] = Vxyz[4][j];eps[4] = Vxyz[5][j];eps[5] = Vxyz[5][j];
+ 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[j] * str[0];
+ Syy += pre[j] * str[1];
+ Szz += pre[j] * str[2];
+ Sxy += pre[j] * str[3];
+ Sxz += pre[j] * str[4];
+ Szy += pre[j] * str[5];
+
+ div += Vxyz[7][j]*dOmega->vpt[j]; /* divergence */
+ vor += Vxyz[8][j]*dOmega->vpt[j]; /* vorticity */
+ }
+
+ }else{
+#endif
for(j=1;j <= vpts;j++) {
/* deviatoric stress, pressure will be added later */
Sxx += 2.0 * pre[j] * (Vxyz[1][j] - dilation[j]); /* */
@@ -371,6 +413,9 @@
div += Vxyz[7][j]*dOmega->vpt[j]; /* divergence */
vor += Vxyz[8][j]*dOmega->vpt[j]; /* vorticity */
}
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ }
+#endif
/* normalize by volume */
Sxx /= E->eco[m][e].area;
Syy /= E->eco[m][e].area;
Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -823,7 +823,12 @@
/* Tracers are placed randomly in cap */
/* (intentionally using rand() instead of srand() )*/
+ /*
+
+ what is this supposed to mean? why not initialize the random
+ number generator with a srand48 call? XXX TWB
+ */
while (E->trace.ntracers[j]<tracers_cap) {
number_of_tries++;
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2010-09-06 05:20:44 UTC (rev 17172)
@@ -48,6 +48,9 @@
static void low_viscosity_channel_factor(struct All_variables *E, float *F);
static void low_viscosity_wedge_factor(struct All_variables *E, float *F);
void parallel_process_termination();
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+#include "anisotropic_viscosity.h"
+#endif
void viscosity_system_input(struct All_variables *E)
@@ -61,6 +64,29 @@
input_boolean("visc_layer_control",&(E->viscosity.layer_control),"off",m);
input_string("visc_layer_file", E->viscosity.layer_file,"visc.dat",m);
+
+ input_boolean("allow_orthotropic_viscosity",&(E->viscosity.allow_orthotropic_viscosity),"off",m);
+#ifndef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ if(E->viscosity.allow_orthotropic_viscosity){ /* error */
+ fprintf(stderr,"error: allow_orthotropic_viscosity is set, but code not compiled with CITCOM_ALLOW_ORTHOTROPIC_VISC\n");
+ parallel_process_termination();
+ }
+#else
+ if(E->viscosity.allow_orthotropic_viscosity){ /* read additional
+ parameters for
+ anisotropic
+ viscosity */
+ input_int("anisotropic_init",&(E->viscosity.anisotropic_init),"0",m);
+ input_string("anisotropic_init_dir",(E->viscosity.anisotropic_init_dir),"",m); /* directory
+ for
+ ggrd
+ type
+ init */
+
+ }
+
+#endif
+ /* allocate memory here */
allocate_visc_vars(E);
for(i=0;i < E->viscosity.num_mat;i++)
@@ -181,41 +207,49 @@
void allocate_visc_vars(struct All_variables *E)
{
- int i;
+ int i,j,k,lim,nel,nno;
- if(E->viscosity.layer_control) {
- /* noz is already defined, but elz is undefined yet */
- i = E->mesh.noz - 1;
- } else {
- i = E->viscosity.num_mat;
- }
+ if(E->viscosity.layer_control) {
+ /* noz is already defined, but elz is undefined yet */
+ i = E->mesh.noz - 1;
+ } else {
+ i = E->viscosity.num_mat;
+ }
+
+ E->viscosity.N0 = (float*) malloc(i*sizeof(float));
+ E->viscosity.E = (float*) malloc(i*sizeof(float));
+ E->viscosity.T = (float*) malloc(i*sizeof(float));
+ E->viscosity.Z = (float*) malloc(i*sizeof(float));
+ E->viscosity.pdepv_a = (float*) malloc(i*sizeof(float));
+ E->viscosity.pdepv_b = (float*) malloc(i*sizeof(float));
+ E->viscosity.pdepv_y = (float*) malloc(i*sizeof(float));
+ E->viscosity.sdepv_expt = (float*) malloc(i*sizeof(float));
+
+ if(E->viscosity.N0 == NULL ||
+ E->viscosity.E == NULL ||
+ E->viscosity.T == NULL ||
+ E->viscosity.Z == NULL ||
+ E->viscosity.pdepv_a == NULL ||
+ E->viscosity.pdepv_b == NULL ||
+ E->viscosity.pdepv_y == NULL ||
+ E->viscosity.sdepv_expt == NULL) {
+ fprintf(stderr, "Error: Cannot allocate visc memory, rank=%d\n",
+ E->parallel.me);
+ parallel_process_termination();
+ }
- E->viscosity.N0 = (float*) malloc(i*sizeof(float));
- E->viscosity.E = (float*) malloc(i*sizeof(float));
- E->viscosity.T = (float*) malloc(i*sizeof(float));
- E->viscosity.Z = (float*) malloc(i*sizeof(float));
- E->viscosity.pdepv_a = (float*) malloc(i*sizeof(float));
- E->viscosity.pdepv_b = (float*) malloc(i*sizeof(float));
- E->viscosity.pdepv_y = (float*) malloc(i*sizeof(float));
- E->viscosity.sdepv_expt = (float*) malloc(i*sizeof(float));
-
- if(E->viscosity.N0 == NULL ||
- E->viscosity.E == NULL ||
- E->viscosity.T == NULL ||
- E->viscosity.Z == NULL ||
- E->viscosity.pdepv_a == NULL ||
- E->viscosity.pdepv_b == NULL ||
- E->viscosity.pdepv_y == NULL ||
- E->viscosity.sdepv_expt == NULL) {
- fprintf(stderr, "Error: Cannot allocate visc memory, rank=%d\n",
- E->parallel.me);
- exit(8);
- }
}
/* ============================================ */
+/*
+ when called from Drive_solvers:
+
+ evisc = E->EVI[E->mesh.levmax]
+ visc = E->VI[E->mesh.levmax]
+
+ */
void get_system_viscosity(E,propogate,evisc,visc)
struct All_variables *E;
int propogate;
@@ -238,7 +272,16 @@
double *TG;
const int vpts = vpoints[E->mesh.nsd];
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ if(E->viscosity.allow_orthotropic_viscosity){
+ if(!E->viscosity.orthotropic_viscosity_init)
+ 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,evisc,propogate);
else
@@ -291,13 +334,32 @@
}
fflush(E->fp_out);
}
-
/* interpolate from gauss quadrature points to node points for output */
visc_from_gint_to_nodes(E,evisc,visc,E->mesh.levmax);
+
if(E->viscosity.SMOOTH){ /* go the other way, for
smoothing */
visc_from_nodes_to_gint(E,visc,evisc,E->mesh.levmax);
}
+
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC /* allow for anisotropy */
+ if(E->viscosity.allow_orthotropic_viscosity){
+ 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);
+
+ if(E->viscosity.SMOOTH){
+ visc_from_nodes_to_gint(E,E->VI2[E->mesh.levmax],E->EVI2[E->mesh.levmax],E->mesh.levmax);
+ visc_from_nodes_to_gint(E,E->VIn1[E->mesh.levmax],E->EVIn1[E->mesh.levmax],E->mesh.levmax);
+ visc_from_nodes_to_gint(E,E->VIn2[E->mesh.levmax],E->EVIn2[E->mesh.levmax],E->mesh.levmax);
+ visc_from_nodes_to_gint(E,E->VIn3[E->mesh.levmax],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
return;
}
Added: mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h (rev 0)
+++ mc/3D/CitcomS/trunk/lib/anisotropic_viscosity.h 2010-09-06 05:20:44 UTC (rev 17172)
@@ -0,0 +1,20 @@
+#ifndef __CITCOM_READ_ANIVISC_HEADER__
+
+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 get_delta(double [3][3][3][3],double [3]);
+void rot_4x4(double [3][3][3][3],double [3][3], double [3][3][3][3]);
+void print_6x6_mat(FILE *, double [6][6]);
+#define __CITCOM_READ_ANIVISC_HEADER__
+#endif
Modified: mc/3D/CitcomS/trunk/lib/ggrd_handling.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/ggrd_handling.h 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/ggrd_handling.h 2010-09-06 05:20:44 UTC (rev 17172)
@@ -45,5 +45,6 @@
void xyz2rtpd(float ,float ,float ,double *);
float find_age_in_MY();
void ggrd_adjust_tbl_rayleigh(struct All_variables *,double **);
+void ggrd_read_anivisc_from_file(struct All_variables *, int );
#define GGRD_DENS_MIN 3200 /* minimum density for PREM scaling of input files */
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2010-09-06 05:20:44 UTC (rev 17172)
@@ -525,7 +525,7 @@
#ifdef USE_GGRD
struct ggrd_master ggrd;
float *surface_rayleigh;
- int ggrd_allow_mixed_vbcs;
+ int ggrd_allow_mixed_vbcs,ggrd_comp_smooth;
float ggrd_vtop_omega[4];
char ggrd_mat_depth_file[1000];
ggrd_boolean ggrd_mat_is_3d;
@@ -783,6 +783,14 @@
float *Vi[NCS],*EVi[NCS];
float *VI[MAX_LEVELS][NCS],*EVI[MAX_LEVELS][NCS];
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ float *VI2[MAX_LEVELS][NCS],*EVI2[MAX_LEVELS][NCS];
+ float *VIn1[MAX_LEVELS][NCS],*EVIn1[MAX_LEVELS][NCS];
+ float *VIn2[MAX_LEVELS][NCS],*EVIn2[MAX_LEVELS][NCS];
+ float *VIn3[MAX_LEVELS][NCS],*EVIn3[MAX_LEVELS][NCS];
+#endif
+
+
int num_zero_resid[MAX_LEVELS][NCS];
int *zero_resid[MAX_LEVELS][NCS];
int *surf_element[NCS],*surf_node[NCS];
Modified: mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2010-09-03 00:55:39 UTC (rev 17171)
+++ mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2010-09-06 05:20:44 UTC (rev 17172)
@@ -35,6 +35,11 @@
int SMOOTH;
int smooth_cycles;
+ int allow_orthotropic_viscosity,orthotropic_viscosity_init;
+#ifdef CITCOM_ALLOW_ORTHOTROPIC_VISC
+ int anisotropic_init; /* 0: isotropic, 1: random, 2: from file */
+ char anisotropic_init_dir[1000];
+#endif
char STRUCTURE[20]; /* which option to determine viscosity field, one of .... */
int FROM_SYSTEM;
More information about the CIG-COMMITS
mailing list