[cig-commits] r5745 - in mc/3D/CitcomS/branches/compressible:
CitcomS/Solver lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Tue Jan 9 14:20:09 PST 2007
Author: tan2
Date: 2007-01-09 14:20:08 -0800 (Tue, 09 Jan 2007)
New Revision: 5745
Modified:
mc/3D/CitcomS/branches/compressible/CitcomS/Solver/Solver.py
mc/3D/CitcomS/branches/compressible/lib/Construct_arrays.c
mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c
mc/3D/CitcomS/branches/compressible/lib/Instructions.c
mc/3D/CitcomS/branches/compressible/lib/Output_h5.c
mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c
mc/3D/CitcomS/branches/compressible/lib/global_defs.h
mc/3D/CitcomS/branches/compressible/module/setProperties.c
Log:
1. Added two non-dimensional number: dissipation # and gruneisen parameter
2. Modified stiffness matrix for TALA
3. Modified stress calculation for TALA
Modified: mc/3D/CitcomS/branches/compressible/CitcomS/Solver/Solver.py
===================================================================
--- mc/3D/CitcomS/branches/compressible/CitcomS/Solver/Solver.py 2007-01-09 22:01:20 UTC (rev 5744)
+++ mc/3D/CitcomS/branches/compressible/CitcomS/Solver/Solver.py 2007-01-09 22:20:08 UTC (rev 5745)
@@ -103,7 +103,7 @@
stream = open("pid%d.cfg" % pid, "w")
else:
stream = None
-
+
self.setProperties(stream)
self.restart = self.inventory.ic.inventory.restart
@@ -282,6 +282,8 @@
datadir_old = inv.str("datadir_old", default="")
rayleigh = inv.float("rayleigh", default=1e+05)
+ dissipation = inv.float("dissipation", default=0.0)
+ gruneisen = inv.float("gruneisen", default=0.0)
Q0 = inv.float("Q0", default=0.0)
stokes_flow_only = inv.bool("stokes_flow_only", default=False)
Modified: mc/3D/CitcomS/branches/compressible/lib/Construct_arrays.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Construct_arrays.c 2007-01-09 22:01:20 UTC (rev 5744)
+++ mc/3D/CitcomS/branches/compressible/lib/Construct_arrays.c 2007-01-09 22:20:08 UTC (rev 5745)
@@ -1,6 +1,6 @@
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
+ *
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,7 +22,7 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
- *
+ *
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#include <math.h>
@@ -341,7 +341,7 @@
for(element=1;element<=nel;element++) {
- get_elt_k(E,element,elt_K,level,m);
+ get_elt_k(E,element,elt_K,level,m,0);
if (E->control.augmented_Lagr)
get_aug_k(E,element,elt_K,level,m);
@@ -633,7 +633,7 @@
for(el=1;el<=E->lmesh.NEL[lev];el++) {
- get_elt_k(E,el,E->elt_k[lev][m][el].k,lev,m); /* not for penalty */
+ get_elt_k(E,el,E->elt_k[lev][m][el].k,lev,m,0);
if (E->control.augmented_Lagr)
get_aug_k(E,el,E->elt_k[lev][m][el].k,lev,m);
Modified: mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c 2007-01-09 22:01:20 UTC (rev 5744)
+++ mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c 2007-01-09 22:20:08 UTC (rev 5745)
@@ -152,21 +152,90 @@
return;
}
+
/*==============================================================
+ Function to supply the element strain-displacement matrix Ba,
+ which is used to compute element stiffness matrix
+ ============================================================== */
+
+void get_ba(struct Shape_function *N, struct Shape_function_dx *GNx,
+ struct CC *cc, struct CCX *ccx, double rtf[4][9],
+ int dims, double ba[9][9][4][7])
+{
+ int k, a, n;
+ const int vpts = VPOINTS3D;
+ const int ends = ENODES3D;
+
+ double ra[9], si[9], ct[9];
+
+ for(k=1;k<=vpts;k++) {
+ ra[k] = rtf[3][k];
+ si[k] = sin(rtf[1][k]);
+ ct[k] = cos(rtf[1][k])/si[k];
+ }
+
+ for(n=1;n<=dims;n++)
+ for(k=1;k<=vpts;k++)
+ for(a=1;a<=ends;a++) {
+ ba[a][k][n][1]=
+ (GNx->vpt[GNVXINDEX(0,a,k)]*cc->vpt[BVINDEX(1,n,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*ccx->vpt[BVXINDEX(1,n,1,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*cc->vpt[BVINDEX(3,n,a,k)]
+ )*ra[k];
+
+ ba[a][k][n][2]=
+ (N->vpt[GNVINDEX(a,k)]*cc->vpt[BVINDEX(1,n,a,k)]*ct[k]
+ + N->vpt[GNVINDEX(a,k)]*cc->vpt[BVINDEX(3,n,a,k)]
+ +(GNx->vpt[GNVXINDEX(1,a,k)]*cc->vpt[BVINDEX(2,n,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*ccx->vpt[BVXINDEX(2,n,2,a,k)]
+ )/si[k]
+ )*ra[k];
+
+ ba[a][k][n][3]=
+ GNx->vpt[GNVXINDEX(2,a,k)]*cc->vpt[BVINDEX(3,n,a,k)];
+
+ ba[a][k][n][4]=
+ (GNx->vpt[GNVXINDEX(0,a,k)]*cc->vpt[BVINDEX(2,n,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*ccx->vpt[BVXINDEX(2,n,1,a,k)]
+ - N->vpt[GNVINDEX(a,k)]*cc->vpt[BVINDEX(2,n,a,k)]*ct[k]
+ +(GNx->vpt[GNVXINDEX(1,a,k)]*cc->vpt[BVINDEX(1,n,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*ccx->vpt[BVXINDEX(1,n,2,a,k)]
+ )/si[k]
+ )*ra[k];
+
+ ba[a][k][n][5]=
+ GNx->vpt[GNVXINDEX(2,a,k)]*cc->vpt[BVINDEX(1,n,a,k)]
+ +(GNx->vpt[GNVXINDEX(0,a,k)]*cc->vpt[BVINDEX(3,n,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*(ccx->vpt[BVXINDEX(3,n,1,a,k)]
+ - cc->vpt[BVINDEX(1,n,a,k)])
+ )*ra[k];
+
+ ba[a][k][n][6]=
+ GNx->vpt[GNVXINDEX(2,a,k)]*cc->vpt[BVINDEX(2,n,a,k)]
+ -ra[k]*N->vpt[GNVINDEX(a,k)]*cc->vpt[BVINDEX(2,n,a,k)]
+ +(GNx->vpt[GNVXINDEX(1,a,k)]*cc->vpt[BVINDEX(3,n,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*ccx->vpt[BVXINDEX(3,n,2,a,k)]
+ )/si[k]*ra[k];
+ }
+
+ return;
+}
+
+/*==============================================================
Function to supply the element k matrix for a given element e.
============================================================== */
-void get_elt_k(E,el,elt_k,lev,m)
+void get_elt_k(E,el,elt_k,lev,m,iconv)
struct All_variables *E;
int el,m;
double elt_k[24*24];
- int lev;
+ int lev, iconv;
{
double bdbmu[4][4];
int pn,qn,ad,bd;
- int a,b,i,j,j1,k,n;
- double rtf[4][9],W[9],ra[9],si[9],ct[9];
+ int a,b,i,j,i1,j1,k;
+ double rtf[4][9],W[9];
void get_global_shape_fn();
void construct_c3x3matrix_el();
struct Shape_function GN;
@@ -174,16 +243,34 @@
struct Shape_function_dx GNx;
double ba[9][9][4][7]; /* integration points,node,3x6 matrix */
+ /* d1 and d2 are the elastic coefficient matrix for incompressible
+ and compressible creeping flow, respectively */
+ const double d1[7][7] =
+ { {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+ {0.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+ {0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 0.0},
+ {0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0} };
+ const double d2[7][7] =
+ { {0., 0., 0., 0., 0., 0., 0.},
+ {0., 2.-2./3., -2./3., -2./3., 0., 0., 0.},
+ {0., -2./3., 2.-2./3., -2./3., 0., 0., 0.},
+ {0., -2./3., -2./3., 2.-2./3., 0., 0., 0.},
+ {0., 0., 0., 0., 1., 0., 0.},
+ {0., 0., 0., 0., 0., 1., 0.},
+ {0., 0., 0., 0., 0., 0., 1.} };
const int nn=loc_mat_size[E->mesh.nsd];
- const int vpts=vpoints[E->mesh.nsd];
- const int ends=enodes[E->mesh.nsd];
+ const int vpts = VPOINTS3D;
+ const int ends = ENODES3D;
const int dims=E->mesh.nsd;
const int sphere_key = 1;
get_global_shape_fn(E,el,&GN,&GNx,&dOmega,0,sphere_key,rtf,lev,m);
- if ((el-1)%E->lmesh.ELZ[lev]==0)
+ 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
@@ -191,70 +278,39 @@
for(k=1;k<=vpts;k++) {
W[k]=g_point[k].weight[dims-1]*dOmega.vpt[k]*E->EVI[lev][m][(el-1)*vpts+k];
+ }
- ra[k] = rtf[3][k];
- si[k] = sin(rtf[1][k]);
- ct[k] = cos(rtf[1][k])/si[k];
- }
+ get_ba(&(E->N), &GNx, &E->element_Cc, &E->element_Ccx,
+ rtf, E->mesh.nsd, ba);
-
- for(k=1;k<=VPOINTS3D;k++)
- for(a=1;a<=ends;a++)
- for(n=1;n<=dims;n++) {
- ba[a][k][n][1]=
- (GNx.vpt[GNVXINDEX(0,a,k)]*E->element_Cc.vpt[BVINDEX(1,n,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Ccx.vpt[BVXINDEX(1,n,1,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Cc.vpt[BVINDEX(3,n,a,k)])*ra[k];
-
- ba[a][k][n][2]=
- (E->N.vpt[GNVINDEX(a,k)]*E->element_Cc.vpt[BVINDEX(1,n,a,k)]*ct[k]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Cc.vpt[BVINDEX(3,n,a,k)]
- +(GNx.vpt[GNVXINDEX(1,a,k)]*E->element_Cc.vpt[BVINDEX(2,n,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Ccx.vpt[BVXINDEX(2,n,2,a,k)])
- /si[k])*ra[k];
-
- ba[a][k][n][3]=
- GNx.vpt[GNVXINDEX(2,a,k)]*E->element_Cc.vpt[BVINDEX(3,n,a,k)];
-
- ba[a][k][n][4]=
- (GNx.vpt[GNVXINDEX(0,a,k)]*E->element_Cc.vpt[BVINDEX(2,n,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Ccx.vpt[BVXINDEX(2,n,1,a,k)]
- - E->N.vpt[GNVINDEX(a,k)]*E->element_Cc.vpt[BVINDEX(2,n,a,k)]*ct[k]
- +(GNx.vpt[GNVXINDEX(1,a,k)]*E->element_Cc.vpt[BVINDEX(1,n,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Ccx.vpt[BVXINDEX(1,n,2,a,k)])
- /si[k])*ra[k];
-
- ba[a][k][n][5]=
- GNx.vpt[GNVXINDEX(2,a,k)]*E->element_Cc.vpt[BVINDEX(1,n,a,k)]
- +(GNx.vpt[GNVXINDEX(0,a,k)]*E->element_Cc.vpt[BVINDEX(3,n,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*(E->element_Ccx.vpt[BVXINDEX(3,n,1,a,k)]
- - E->element_Cc.vpt[BVINDEX(1,n,a,k)]))*ra[k];
-
- ba[a][k][n][6]=
- GNx.vpt[GNVXINDEX(2,a,k)]*E->element_Cc.vpt[BVINDEX(2,n,a,k)]
- -ra[k]*E->N.vpt[GNVINDEX(a,k)]*E->element_Cc.vpt[BVINDEX(2,n,a,k)]
- +(GNx.vpt[GNVXINDEX(1,a,k)]*E->element_Cc.vpt[BVINDEX(3,n,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Ccx.vpt[BVXINDEX(3,n,2,a,k)])
- /si[k]*ra[k];
- }
-
for(a=1;a<=ends;a++)
for(b=a;b<=ends;b++) {
- ad=dims*(a-1);
- bd=dims*(b-1);
-
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.gruneisen > 0)
for(i=1;i<=dims;i++)
for(j=1;j<=dims;j++)
- for(j1=1;j1<=6;j1++)
- for(k=1;k<=VPOINTS3D;k++)
- bdbmu[i][j] +=
- W[k]*((j1>3)?1.0:2.0)*ba[a][k][i][j1]*ba[b][k][j][j1];
+ for(i1=1;i1<=6;i1++)
+ for(j1=1;j1<=6;j1++)
+ for(k=1;k<=VPOINTS3D;k++)
+ bdbmu[i][j] +=
+ W[k]*d2[i1][j1]*ba[a][k][i][i1]*ba[b][k][j][j1];
+ else
+ for(i=1;i<=dims;i++)
+ for(j=1;j<=dims;j++)
+ for(i1=1;i1<=6;i1++)
+ for(j1=1;j1<=6;j1++)
+ for(k=1;k<=VPOINTS3D;k++)
+ bdbmu[i][j] +=
+ W[k]*d1[i1][j1]*ba[a][k][i][i1]*ba[b][k][j][j1];
+
/**/
+ ad=dims*(a-1);
+ bd=dims*(b-1);
+
pn=ad*nn+bd;
qn=bd*nn+ad;
@@ -793,7 +849,7 @@
nodeb=E->ien[m][el].node[b];
if ((E->node[m][nodeb]&type)&&(E->sphere.cap[m].VB[j][nodeb]!=0.0)){
if(!got_elt_k) {
- get_elt_k(E,el,elt_k,E->mesh.levmax,m);
+ get_elt_k(E,el,elt_k,E->mesh.levmax,m,1);
got_elt_k = 1;
}
q = dims*(b-1)+j-1;
Modified: mc/3D/CitcomS/branches/compressible/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Instructions.c 2007-01-09 22:01:20 UTC (rev 5744)
+++ mc/3D/CitcomS/branches/compressible/lib/Instructions.c 2007-01-09 22:20:08 UTC (rev 5745)
@@ -364,6 +364,8 @@
input_int("piterations",&(E->control.p_iterations),"100,0,nomax",m);
input_float("rayleigh",&(E->control.Atemp),"essential",m);
+ input_float("dissipation",&(E->control.Di),"0.0",m);
+ input_float("gruneisen",&(E->control.gruneisen),"0.0",m);
/* data section */
input_float("Q0",&(E->control.Q0),"0.0",m);
Modified: mc/3D/CitcomS/branches/compressible/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Output_h5.c 2007-01-09 22:01:20 UTC (rev 5744)
+++ mc/3D/CitcomS/branches/compressible/lib/Output_h5.c 2007-01-09 22:20:08 UTC (rev 5745)
@@ -1458,6 +1458,8 @@
status = set_attribute_string(input, "datafile_old", E->control.old_P_file);
status = set_attribute_float(input, "rayleigh", E->control.Atemp);
+ status = set_attribute_float(input, "dissipation", E->control.Di);
+ status = set_attribute_float(input, "gruneisen", E->control.gruneisen);
status = set_attribute_float(input, "Q0", E->control.Q0);
status = set_attribute_int(input, "stokes_flow_only", E->control.stokes);
Modified: mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c 2007-01-09 22:01:20 UTC (rev 5744)
+++ mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c 2007-01-09 22:20:08 UTC (rev 5745)
@@ -191,6 +191,7 @@
int i,j,e,node,m;
float VV[4][9],Vxyz[9][9],Szz,Sxx,Syy,Sxy,Sxz,Szy,div,vor;
+ double dilation[9];
double pre[9],tww[9],rtf[4][9];
double velo_scaling, stress_scaling;
@@ -220,6 +221,7 @@
for(j=1;j<=vpts;j++) {
pre[j] = E->EVi[m][(e-1)*vpts+j]*dOmega.vpt[j];
+ dilation[j] = 0.0;
Vxyz[1][j] = 0.0;
Vxyz[2][j] = 0.0;
Vxyz[3][j] = 0.0;
@@ -262,16 +264,24 @@
+ VV[2][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j])
- VV[1][i]*GNx.vpt[GNVXINDEX(1,i,j)] );
}
- Sxx += 2.0 * pre[j] * Vxyz[1][j];
- Syy += 2.0 * pre[j] * Vxyz[2][j];
- Szz += 2.0 * pre[j] * Vxyz[3][j];
- Sxy += pre[j] * Vxyz[4][j];
- Sxz += pre[j] * Vxyz[5][j];
- Szy += pre[j] * Vxyz[6][j];
- div += Vxyz[7][j]*dOmega.vpt[j];
- vor += Vxyz[8][j]*dOmega.vpt[j];
}
+ if(E->control.gruneisen > 0) {
+ for(j=1;j<=vpts;j++)
+ dilation[j] = (Vxyz[1][j] + Vxyz[2][j] + Vxyz[3][j]) * 2.0 / 3.0;
+ }
+
+ for(j=1;j<=vpts;j++) {
+ Sxx += 2.0 * pre[j] * (Vxyz[1][j] - dilation[j]);
+ Syy += 2.0 * pre[j] * (Vxyz[2][j] - dilation[j]);
+ Szz += 2.0 * pre[j] * (Vxyz[3][j] - dilation[j]);
+ Sxy += pre[j] * Vxyz[4][j];
+ Sxz += pre[j] * Vxyz[5][j];
+ Szy += pre[j] * Vxyz[6][j];
+ div += Vxyz[7][j]*dOmega.vpt[j];
+ vor += Vxyz[8][j]*dOmega.vpt[j];
+ }
+
Sxx /= E->eco[m][e].area;
Syy /= E->eco[m][e].area;
Szz /= E->eco[m][e].area;
@@ -676,8 +686,8 @@
/* eu [m*dims+2] = VV[3][m+1]; */
/* } */
-/* get_elt_k(E,el,eltk,lev,j); */
-/* get_elt_k(E,elb,eltkb,lev,j); */
+/* get_elt_k(E,el,eltk,lev,j,1); */
+/* get_elt_k(E,elb,eltkb,lev,j,1); */
/* get_elt_f(E,el,eltf,0,j); */
/* get_elt_f(E,elb,eltfb,0,j); */
/* get_elt_g(E,el,eltg,lev,j); */
Modified: mc/3D/CitcomS/branches/compressible/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/global_defs.h 2007-01-09 22:01:20 UTC (rev 5744)
+++ mc/3D/CitcomS/branches/compressible/lib/global_defs.h 2007-01-09 22:20:08 UTC (rev 5745)
@@ -468,6 +468,7 @@
int post_p;
int post_topo;
int SLAB;
+ int compressible;
char GEOMETRY[20]; /* one of ... */
int CART2D;
@@ -515,6 +516,9 @@
float sob_tolerance;
float Atemp;
+ float Di;
+ float gruneisen;
+
float inputdiff;
float VBXtopval;
float VBXbotval;
Modified: mc/3D/CitcomS/branches/compressible/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/module/setProperties.c 2007-01-09 22:01:20 UTC (rev 5744)
+++ mc/3D/CitcomS/branches/compressible/module/setProperties.c 2007-01-09 22:20:08 UTC (rev 5745)
@@ -450,6 +450,8 @@
getStringProperty(properties, "datafile_old", E->control.data_prefix_old, fp);
getFloatProperty(properties, "rayleigh", E->control.Atemp, fp);
+ getFloatProperty(properties, "dissipation", E->control.Di, fp);
+ getFloatProperty(properties, "gruneisen", E->control.gruneisen, fp);
getFloatProperty(properties, "Q0", E->control.Q0, fp);
getIntProperty(properties, "stokes_flow_only", E->control.stokes, fp);
More information about the cig-commits
mailing list