[cig-commits] r17256 - in mc/3D/CitcomCU/trunk/src: . expokit
becker at geodynamics.org
becker at geodynamics.org
Sun Oct 10 21:21:33 PDT 2010
Author: becker
Date: 2010-10-10 21:21:32 -0700 (Sun, 10 Oct 2010)
New Revision: 17256
Added:
mc/3D/CitcomCU/trunk/src/expokit/
mc/3D/CitcomCU/trunk/src/expokit/dgpadm.f
Modified:
mc/3D/CitcomCU/trunk/src/Advection_diffusion.c
mc/3D/CitcomCU/trunk/src/Anisotropic_viscosity.c
mc/3D/CitcomCU/trunk/src/Boundary_conditions.c
mc/3D/CitcomCU/trunk/src/Citcom.c
mc/3D/CitcomCU/trunk/src/Composition_adv.c
mc/3D/CitcomCU/trunk/src/Construct_arrays.c
mc/3D/CitcomCU/trunk/src/Drive_solvers.c
mc/3D/CitcomCU/trunk/src/Element_calculations.c
mc/3D/CitcomCU/trunk/src/General_matrix_functions.c
mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
mc/3D/CitcomCU/trunk/src/Global_operations.c
mc/3D/CitcomCU/trunk/src/Instructions.c
mc/3D/CitcomCU/trunk/src/Makefile
mc/3D/CitcomCU/trunk/src/Makefile.gzdir
mc/3D/CitcomCU/trunk/src/Makefile.gzdir_ani
mc/3D/CitcomCU/trunk/src/Nodal_mesh.c
mc/3D/CitcomCU/trunk/src/Output_gzdir.c
mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c
mc/3D/CitcomCU/trunk/src/Size_does_matter.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/anisotropic_viscosity.h
mc/3D/CitcomCU/trunk/src/global_defs.h
mc/3D/CitcomCU/trunk/src/prototypes.h
mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h
Log:
Further modified anisotropic viscosity. Program stalls on our cluster
for some combinations of larger (>32) number of CPUs for reasons
unknown, and unclear how this related to any recent modifications.
Modified: mc/3D/CitcomCU/trunk/src/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Advection_diffusion.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Advection_diffusion.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -141,7 +141,6 @@
{
}
-
if(on_off == 0)
{
E->advection.timesteps++;
@@ -187,7 +186,6 @@
corrector(E, E->T, E->Tdot, DTdot);
}
}
-
/* get the max temperature for new T */
E->monitor.T_interior = Tmax(E, E->T);
@@ -209,18 +207,15 @@
count++;
temperatures_conform_bcs(E);
-
E->advection.last_sub_iterations = count;
Euler(E, E->C, E->V, on_off);
-
+
E->monitor.elapsed_time += E->advection.timestep;
} /* end for on_off==0 */
-
thermal_buoyancy(E);
-
if(E->monitor.solution_cycles < E->advection.max_timesteps)
E->control.keep_going = 1;
else
Modified: mc/3D/CitcomCU/trunk/src/Anisotropic_viscosity.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Anisotropic_viscosity.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Anisotropic_viscosity.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -1,5 +1,7 @@
/*
+ set of routines to deal with anisotropic viscosity
+
orthotropic viscosity following Muehlhaus, Moresi, Hobbs and Dufour
(PAGEOPH, 159, 2311, 2002)
@@ -28,20 +30,33 @@
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 */
+ (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 */
+ }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);
+ if(E->avmode[lev][off] == CITCOM_ANIVISC_ORTHO_MODE){
+ /*
+ orthotropic
+ */
+ get_constitutive_orthotropic_viscosity(D,E->EVI2[lev][off],
+ n,convert_to_spherical,theta,phi);
+ }else if(E->avmode[lev][off] == CITCOM_ANIVISC_TI_MODE){
+ /*
+ transversely isotropic
+ */
+ get_constitutive_ti_viscosity(D,E->EVI2[lev][off],0.,
+ n,convert_to_spherical,
+ theta,phi);
}
+ //if(E->EVI2[lev][off] != 0) print_6x6_mat(stderr,D);
}
}else{
get_constitutive_isotropic(D);
@@ -218,81 +233,105 @@
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;
- }
+ if(init_stage){
+ if(E->parallel.me == 0)
+ fprintf(stderr,"set_anisotropic_viscosity: allowing for %s viscosity\n",
+ (E->viscosity.allow_anisotropic_viscosity == 1)?("orthotropic"):("transversely isotropic"));
+ 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 3: /* first init for vel */
+ case 4: /* ISA */
+ case 5: /* same for mixed alignment */
+ 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;
+ E->avmode[i][off] = (unsigned char)
+ E->viscosity.allow_anisotropic_viscosity;
}
}
- break;
- case 1: /*
+ }
+ 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];
- }
+ 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++){
+ vis2 = 0.01+drand48()*0.99; /* random fluctuation */
+ /* 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];
+ E->avmode[i][off] = (unsigned char)
+ E->viscosity.allow_anisotropic_viscosity;
+
}
}
- break;
- case 2: /* from file */
+ }
+ 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();
+ 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);
+ 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{
+ //if(E->parallel.me == 0)fprintf(stderr,"reassigning anisotropic viscosity, mode %i\n",E->viscosity.anisotropic_init);
+ if((E->monitor.solution_cycles > 0) || (E->monitor.visc_iter_count > 0)){
+ /* standard operation every time step */
+ switch(E->viscosity.anisotropic_init){
+ /* if flow has been computed, use director aligned with ISA */
+ case 3: /* vel lignment */
+ /* we have a velocity solution already */
+
+ align_director_with_ISA_for_element(E,CITCOM_ANIVISC_ALIGN_WITH_VEL);
break;
- default:
- fprintf(stderr,"set_anisotropic_viscosity_at_element_level: anisotropic_init %i undefined\n",
- E->viscosity.anisotropic_init);
- parallel_process_termination();
+ case 4: /* vel alignment */
+ if((E->monitor.solution_cycles > 0) || (E->monitor.visc_iter_count > 0))
+ align_director_with_ISA_for_element(E,CITCOM_ANIVISC_ALIGN_WITH_ISA);
break;
+ case 5: /* mixed alignment */
+ if((E->monitor.solution_cycles > 0) || (E->monitor.visc_iter_count > 0))
+ align_director_with_ISA_for_element(E,CITCOM_ANIVISC_MIXED_ALIGN);
+ break;
+ default: /* default, no action */
+ break;
}
- E->viscosity.anisotropic_viscosity_init = TRUE;
- /* end initialization stage */
- }else{
- /* standard operation every time step */
-
-
}
- } /* end anisotropic viscosity branch */
+ }
}
@@ -327,7 +366,7 @@
*/
- void conv_cart4x4_to_spherical(double c[3][3][3][3], double theta, double phi, double p[3][3][3][3])
+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);
@@ -359,7 +398,7 @@
*/
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],
+ double a[3][3][3][3],b[3][3][3][3],rot[3][3],
hlen2,x2,y2,xy,zm1;
/* normalize */
normalize_vec3d((n+0),(n+1),(n+2));
@@ -411,7 +450,8 @@
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) */
+ 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);
@@ -423,13 +463,209 @@
}
+/*
+ mode =
+ CITCOM_ANIVISC_ALIGN_WITH_VEL: align with velocity
+ CITCOM_ANIVISC_ALIGN_WITH_ISA: align with ISA
+ CITCOM_ANIVISC_MIXED_ALIGN: mixed alignment
+
+ */
+void align_director_with_ISA_for_element(struct All_variables *E,
+ int mode)
+{
+ double rtf[4][9];
+ double VV[4][9],lgrad[3][3],isa[3],evel[3];
+ int e,i,off;
+ float base[9],n[3];
+ static struct CC Cc;
+ static struct CCX Ccx;
+ const int dims = E->mesh.nsd;
+ const int ppts = ppoints[dims];
+ const int vpts = vpoints[dims];
+ const int ends = enodes[dims];
+ const int lev = E->mesh.levmax;
+ const int nel = E->lmesh.nel;
+ unsigned char avmode;
+ double vis2 ;
+ /* anisotropy maximum strength */
+ vis2 = 1. - E->viscosity.ani_vis2_factor; /* 1-eta_w/eta_s */
+
+ if(E->parallel.me == 0){
+ switch(mode){
+ case CITCOM_ANIVISC_ALIGN_WITH_VEL:
+ fprintf(stderr,"align_director_with_ISA_for_element: aligning, max ani %g, align with vel\n",
+ vis2);
+ break;
+ case CITCOM_ANIVISC_ALIGN_WITH_ISA:
+ fprintf(stderr,"align_director_with_ISA_for_element: aligning, max ani %g, align with ISA\n",
+ vis2);
+ break;
+ case CITCOM_ANIVISC_MIXED_ALIGN:
+ fprintf(stderr,"align_director_with_ISA_for_element: aligning, max ani %g, align with uniaxial/vel\n",
+ vis2);
+ break;
+ default:
+ fprintf(stderr,"align_director_with_ISA_for_element: mode %i undefined\n",mode);
+ myerror("",E);
+ }
+ }
+ for(e=1; e <= nel; e++) {
+ if(((E->viscosity.anivisc_layer > 0)&&
+ (E->mat[e] <= E->viscosity.anivisc_layer))||
+ ((E->viscosity.anivisc_layer < 0)&&
+ (E->mat[e] == -E->viscosity.anivisc_layer))){
+
+ if(E->control.Rsphere){ /* need rtf for spherical */
+ get_rtf(E, e, 1, rtf, lev); /* pressure points */
+ //if((e-1)%E->lmesh.elz==0)
+ construct_c3x3matrix_el(E,e,&Cc,&Ccx,lev,1);
+ }
+ for(i = 1; i <= ends; i++){ /* velocity at element nodes */
+ off = E->ien[e].node[i];
+ VV[1][i] = E->V[1][off];
+ VV[2][i] = E->V[2][off];
+ VV[3][i] = E->V[3][off];
+ }
+ /* calculate velocity gradient matrix */
+ get_vgm_p(VV,&(E->N),&(E->GNX[lev][e]),&Cc, &Ccx,rtf,
+ E->mesh.nsd,ppts,ends,(E->control.Rsphere),lgrad,
+ evel);
+ /* calculate the ISA axis and determine the type of
+ anisotropy */
+ avmode = calc_isa_from_vgm(lgrad,evel,e,isa,E,mode);
+ /*
+ convert for spherical (Citcom system) to Cartesian
+ */
+ if(E->control.Rsphere){
+ calc_cbase_at_tp(rtf[1][1],rtf[2][1],base);
+ convert_pvec_to_cvec(isa[2],isa[0],isa[1],base,n);
+ }else{
+ n[0]=isa[0];n[1]=isa[1];n[2]=isa[2];
+ }
+ /* assign to director for all vpoints */
+ for(i=1;i <= vpts;i++){
+ off = (e-1)*vpts + i;
+ E->avmode[lev][off] = avmode;
+ E->EVI2[lev][off] = vis2;
+ E->EVIn1[lev][off] = n[0];
+ E->EVIn2[lev][off] = n[1];
+ E->EVIn3[lev][off] = n[2];
+ }
+ } /* end in layer branch */
+ }
+}
/*
+ compute the ISA axis from velocity gradient tensor l, element
+ velocity evel, and element number e
+
+ mode input:
+ CITCOM_ANIVISC_ALIGN_WITH_VEL: align with velocity
+ CITCOM_ANIVISC_ALIGN_WITH_ISA: align with ISA
+ CITCOM_ANIVISC_MIXED_ALIGN: mixed alignment
+
+ returns type of anisotropy CITCOM_ANIVISC_ORTHO_MODE : orthotropic (shear/normal)
+ CITCOM_ANIVISC_TI_MODE : TI
+
+*/
+unsigned char calc_isa_from_vgm(double l[3][3], double ev[3],
+ int e,double isa[3],
+ struct All_variables *E,
+ int mode)
+{
+ double d[3][3],r[3][3],ltrace,eval[3],lc[3][3],
+ evec[3][3],strain_scale,gol,t1,t2;
+ int i,j,isa_mode;
+ /* copy */
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++)
+ lc[i][j] = l[i][j];
+ remove_trace_3x3(lc);
+ calc_strain_from_vgm(lc,d); /* strain-rate */
+#ifndef USE_GGRD
+ myerror("need USE_GGRD compile for ISA axes",E);
+#else
+ ggrd_solve_eigen3x3(d,eval,evec,E); /* compute the eigensystem */
+#endif
+ /* normalize by largest abs(eigenvalue) */
+ strain_scale = ((t1=fabs(eval[2])) > (t2=fabs(eval[0])))?(t1):(t2);
+ for(i=0;i<3;i++){
+ eval[i] /= strain_scale; /* normalized eigenvalues */
+ for(j=0;j<3;j++)
+ lc[i][j] /= strain_scale;
+ }
+ /* recompute normalized strain rate and rotation */
+ calc_strain_from_vgm(lc,d);
+ calc_rot_from_vgm(lc,r); /* rotation */
+ switch(mode){
+ case CITCOM_ANIVISC_ALIGN_WITH_VEL: /* use velocity */
+ isa[0] = ev[0];isa[1] = ev[1];isa[2]=ev[2];
+ normalize_vec3d(isa,(isa+1),(isa+2));
+ return CITCOM_ANIVISC_TI_MODE;
+ break;
+ case CITCOM_ANIVISC_ALIGN_WITH_ISA:
+ isacalc(lc,&gol,isa,E,&isa_mode);
+ if((isa_mode == -1)||(isa_mode == 0)){
+ /* ISA cannot be found = align with flow */
+ isa[0] = ev[0];isa[1] = ev[1];isa[2]=ev[2];
+ normalize_vec3d(isa,(isa+1),(isa+2));
+ return CITCOM_ANIVISC_TI_MODE;
+ }else{ /* actual ISA */
+ return CITCOM_ANIVISC_ORTHO_MODE;
+ }
+ break;
+ case CITCOM_ANIVISC_MIXED_ALIGN:
+ /* mixed */
+ if(is_pure_shear(lc,d,r)){
+ /* align any pure shear state with the biggest absolute
+ eigenvector */
+ /* largest EV (extensional) */
+ //isa[0] = evec[0][0];isa[1] = evec[0][1];isa[2] = evec[0][2];
+ /* smallest (compressive) EV */
+ isa[0] = evec[2][0];isa[1] = evec[2][1];isa[2] = evec[2][2];
+ return CITCOM_ANIVISC_ORTHO_MODE;
+ }else{
+ /* simple shear */
+ isa[0] = ev[0];isa[1] = ev[1];isa[2]=ev[2]; /* align with vel for now */
+ normalize_vec3d(isa,(isa+1),(isa+2));
+ return CITCOM_ANIVISC_TI_MODE; /* TI */
+ }
+ break;
+ default:
+ myerror("ISA mode undefined",E);
+ break;
+ }
+ return 0;
+}
+/*
+ determine if deformation state is pure shear from normalized
+ velocity gradient, strain, and rotation tensor
+
+*/
+
+int is_pure_shear(double l[3][3],double e[3][3],double r[3][3])
+{
+
+ double mrot,tmp,e2;
+ /* find max rotation */
+ mrot = fabs(r[0][1]); /* xy */
+ if((tmp=fabs(r[0][2]))>mrot)mrot = tmp; /* xz */
+ if((tmp=fabs(r[1][2]))>mrot)mrot = tmp; /* yz */
+ /* second invariant of strain-rate tensor */
+ e2 = second_invariant_from_3x3(e);
+ if(mrot/e2 >= 1)
+ return 0; /* simple shear */
+ else
+ return 1; /* pure shear */
+
+}
+/*
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])
+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;
@@ -444,8 +680,9 @@
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];
+ 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])
@@ -490,10 +727,19 @@
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,"%14.5e ",(fabs(c[i][j])<5e-15)?(0):(c[i][j]));
fprintf(out,"\n");
}
}
+void print_3x3_mat(FILE *out, double c[3][3])
+{
+ int i,j;
+ for(i=0;i<3;i++){
+ for(j=0;j<3;j++)
+ fprintf(out,"%14.5e ",(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
@@ -590,36 +836,201 @@
c[j][i] = c[i][j];
}
-void mat_mult_vec_3x3(double a[3][3],double b[3],double c[3])
+void isacalc(double l[3][3], double *gol,double isa[3],
+ struct All_variables *E,int *isa_mode)
{
- 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];
+
+ //
+ // input: l: *normalized* velocity gradient tensor, in FORTRAN sorting
+ // output: isa(3): infite strain axis,
+ // gol: grain orientation lag
+ double ltrace,isa_diff;
+ double f[3][3],eval[3],evec[3][3],u[3][3],le[3][3];
+
+ // make sure l does not have a trace
+ remove_trace_3x3(l);
+
+#ifdef CITCOM_USE_EXPOKIT
+ calc_exp_matrixt(l,75,le,E); /* following Kaminski & Ribe (G-Cubed, 2001) */
+ f_times_ft(le,u);ggrd_solve_eigen3x3(u,eval,evec,E);
+ isa[0] = evec[0][0]; isa[1] = evec[0][1]; isa[2] = evec[0][2];
+ calc_exp_matrixt(l,80,le,E); f_times_ft(le,u);ggrd_solve_eigen3x3(u,eval,evec,E);
+ isa_diff = 1.-fabs(evec[0][0]*isa[0]+evec[0][1]*isa[1]+evec[0][2]*isa[2]);
+ if(isa_diff > 1e-4) /* ISA does not exist */
+ *isa_mode = -1;
+ else
+ *isa_mode = 1;
+ //fprintf(stderr,"A: %11g %11g %11g - %11g %i\n",isa[0],isa[1],isa[2],isa_diff,*isa_mode);
+
+#else
+ // Limit deformation gradient tensor for infinite time
+ // calculation of the ISE orientation using Sylvester's formula
+ drex_eigen(l,f,isa_mode);
+ if(*isa_mode == -1){
+ // isa is flow1
+ isa[0]=isa[1]=isa[2] = -1.;
+ }else if(*isa_mode == 0){
+ isa[0] = isa[1] = isa[2] = 0;
+ *gol = -1.;
+ }else{
+ // 2. formation of the left-stretch tensor U = FFt
+ f_times_ft(f,u);
+ // 3. eigen-values and eigen-vectors of U
+ ggrd_solve_eigen3x3(u,eval,evec,E);
+ // largest eigenvector
+ isa[0] = evec[0][0];isa[1] = evec[0][1];isa[2] = evec[0][2];
+ //fprintf(stderr,"B: %11g %11g %11g\n",evec[0][0],evec[0][1],evec[0][2]);
}
+#endif
}
-void normalize_vec3(float *x, float *y, float *z)
+/*
+ F^2 = F * F^T
+ */
+void f_times_ft(double f[3][3],double out[3][3])
{
- 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;
+ int i,j,k;
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++){
+ out[i][j] = 0.0;
+ for(k=0;k<3;k++)
+ out[i][j] += f[i][k] * f[j][k];
+ }
}
-void normalize_vec3d(double *x, double *y, double *z)
+/* find eigenvalues of velocity gradient tensor (modified from DREX
+ code of Kaminski et al. 2004)
+
+*/
+void drex_eigen(double l[3][3],double f[3][3], int *mode)
{
- double len = 0.;
- len += (*x) * (*x);len += (*y) * (*y);len += (*z) * (*z);
- len = sqrt(len);
- *x /= len;*y /= len;*z /= len;
+ double a2,a3,q,q3,r2,r,theta,xx,lambda1,lambda2,lambda3,sq2;
+ const double four_pi = 4.0*M_PI,
+ two_pi = 2.0*M_PI,
+ one_third=1./3.;
+ /* looking for the eigen-values of L (using tr(l)=0) */
+ a2 = l[0][0] * l[1][1] + l[1][1] * l[2][2] + l[2][2]*l[0][0] -
+ l[0][1] * l[1][0] - l[1][2]*l[2][1] - l[2][0]*l[0][2];
+ a3 = l[0][0]*l[1][2]*l[2][1] + l[0][1]*l[1][0]*l[2][2] +
+ l[0][2]*l[1][1]*l[2][0] - l[0][0]*l[1][1]*l[2][2] -
+ l[0][1]*l[1][2]*l[2][0] - l[0][2]*l[1][0]*l[2][1];
+
+ q = -a2/3.;
+ r = a3/2.;
+ q3 = q*q*q;
+ r2 = r*r;
+
+ if(fabs(q) < 1e-9){
+ /* simple shear, isa=veloc */
+ *mode = -1;
+ }else if(q3-r2 >= 0){
+ sq2 = 2*sqrt(q);
+
+ theta = acos(pow(r/q,1.5));
+ lambda1 = -sq2*cos(theta/3);
+ lambda2 = -sq2*cos((theta+two_pi)/3.);
+ lambda3 = -sq2*cos((theta+four_pi)/3.);
+
+ if (fabs(lambda1-lambda2) < 1e-13)
+ lambda1 = lambda2;
+ if (fabs(lambda2-lambda3) < 1e-13)
+ lambda2 = lambda3;
+ if (fabs(lambda3-lambda1) < 1e-13)
+ lambda3 = lambda1;
+
+ if((lambda1 > lambda2) && (lambda1 > lambda3)) {
+ malmul_scaled_id(f,l,-lambda2,-lambda3);*mode=1;
+ }else if((lambda2 > lambda3 ) && (lambda2 > lambda1)){
+ malmul_scaled_id(f,l,-lambda3,-lambda1);*mode = 1;
+ }else if((lambda3 > lambda1 ) && (lambda3 > lambda2)){
+ malmul_scaled_id(f,l,-lambda1,-lambda2);*mode = 1;
+ }else if((lambda1 == lambda2 ) && (lambda3 > lambda1)) {
+ malmul_scaled_id(f,l,-lambda1,-lambda2);*mode = 1;
+ }else if((lambda2 == lambda3 ) && (lambda1 > lambda2)){
+ malmul_scaled_id(f,l,-lambda2,-lambda3);*mode = 1;
+ }else if((lambda3 == lambda1 ) && (lambda2 > lambda3)){
+ malmul_scaled_id(f,l,-lambda3,-lambda1);*mode = 1;
+ }else if((lambda1 == lambda2 ) && (lambda3 < lambda1)){
+ *mode =0;
+ }else if((lambda2 == lambda3 ) && (lambda1 < lambda2)){
+ *mode = 0;
+ }else if((lambda3 == lambda1 ) && (lambda2 < lambda3)){
+ *mode = 0;
+ }
+ }else{
+ xx = pow(sqrt(r2-q3)+fabs(r),one_third);
+ if(r < 0)
+ lambda1 = xx+q/xx;
+ else
+ lambda1 = -xx+q/xx;
+ lambda2 = -lambda1/2.;
+ lambda3 = -lambda1/2.;
+ if (lambda1 > 1e-9) {
+ malmul_scaled_id(f,l,-lambda2,-lambda3);
+ *mode = 2;
+ }else{
+ *mode = 0;
+ }
+
+ }
}
-void cross_product(double a[3],double b[3],double c[3])
+/* f = matmul(l-lambda2*Id,l-lambda3*Id); */
+void malmul_scaled_id(double f[3][3],double l[3][3],
+ double f1,double f2)
{
- 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];
+ double a[3][3],b[3][3];
+ int i,j;
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++){
+ a[i][j] = b[i][j] = l[i][j];
+ if(i==j){
+ a[i][j] += f1;
+ b[i][j] += f2;
+ }
+ }
+ matmul_3x3(a,b,f);
+
+
}
+#ifdef CITCOM_USE_EXPOKIT
+/*
+calculate exp(A t) using DGPADM from EXPOKIT for a 3x3
+
+(needs LAPACK)
+
+*/
+void calc_exp_matrixt(double a[3][3],double t,double ae[3][3],
+ struct All_variables *E)
+{
+ int ideg=6;// degre of Pade approximation, six should be ok
+ int m=3,ldh=3;// a(ldh,m) dimensions
+ double *wsp;// work space
+ double af[9];
+ int ipiv[4],iexph,ns;// workspace, output pointer, nr of squareas
+ int iflag,lwsp;// exit code, size of workspace
+ int i,j,k;
+ /*
+ work space
+ */
+ lwsp = 2*(4*m*m+ideg+1);// size of workspace, oversized by factor two
+ wsp = (double *)calloc(lwsp,sizeof(double));
+ /* assign fortran style */
+ for(k=i=0;i<3;i++)
+ for(j=0;j<3;j++,k++)
+ af[k] = a[i][j];
+ //
+ // call to expokit routine
+ dgpadm_(&ideg,&m,&t,af,&ldh,wsp,&lwsp,ipiv,&iexph,&ns,&iflag);
+ if(iflag < 0)
+ myerror("calc_exp_matrixt: problem in dgpadm",E);
+ // assign to output
+ for(i=0,k=iexph-1;i<3;i++)
+ for(j=0;j<3;j++,k++)
+ ae[i][j] = wsp[k];
+
+ free(wsp);
+
+}
#endif
+
+
+#endif
Modified: mc/3D/CitcomCU/trunk/src/Boundary_conditions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Boundary_conditions.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Boundary_conditions.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -44,28 +44,28 @@
void velocity_boundary_conditions(struct All_variables *E)
{
int lv;
- //int node, d;
- int node;
-
+ int node,bottom,top;
+ bottom = 1;
for(lv = E->mesh.levmax; lv >= E->mesh.levmin; lv--)
{
- if(E->mesh.botvbc != 1)
+ top = E->mesh.NOZ[lv];
+ if(E->mesh.botvbc != 1) /* free slip on bottom */
{
- horizontal_bc(E, E->VB, 1, 1, 0.0, VBX, 0, lv);
- horizontal_bc(E, E->VB, 1, 3, 0.0, VBZ, 1, lv);
- horizontal_bc(E, E->VB, 1, 2, 0.0, VBY, 0, lv);
- horizontal_bc(E, E->VB, 1, 1, E->control.VBXbotval, SBX, 1, lv);
- horizontal_bc(E, E->VB, 1, 3, 0.0, SBZ, 0, lv);
- horizontal_bc(E, E->VB, 1, 2, E->control.VBYbotval, SBY, 1, lv);
+ horizontal_bc(E, E->VB, bottom, 1, 0.0, VBX, 0, lv);
+ horizontal_bc(E, E->VB, bottom, 3, 0.0, VBZ, 1, lv);
+ horizontal_bc(E, E->VB, bottom, 2, 0.0, VBY, 0, lv);
+ horizontal_bc(E, E->VB, bottom, 1, E->control.VBXbotval, SBX, 1, lv);
+ horizontal_bc(E, E->VB, bottom, 3, 0.0, SBZ, 0, lv);
+ horizontal_bc(E, E->VB, bottom, 2, E->control.VBYbotval, SBY, 1, lv);
}
- if(E->mesh.topvbc != 1)
+ if(E->mesh.topvbc != 1) /* free slip on top */
{
- horizontal_bc(E, E->VB, E->mesh.NOZ[lv], 1, 0.0, VBX, 0, lv);
- horizontal_bc(E, E->VB, E->mesh.NOZ[lv], 3, 0.0, VBZ, 1, lv);
- horizontal_bc(E, E->VB, E->mesh.NOZ[lv], 2, 0.0, VBY, 0, lv);
- horizontal_bc(E, E->VB, E->mesh.NOZ[lv], 1, E->control.VBXtopval, SBX, 1, lv);
- horizontal_bc(E, E->VB, E->mesh.NOZ[lv], 3, 0.0, SBZ, 0, lv);
- horizontal_bc(E, E->VB, E->mesh.NOZ[lv], 2, E->control.VBYtopval, SBY, 1, lv);
+ horizontal_bc(E, E->VB, top, 1, 0.0, VBX, 0, lv);
+ horizontal_bc(E, E->VB, top, 3, 0.0, VBZ, 1, lv);
+ horizontal_bc(E, E->VB, top, 2, 0.0, VBY, 0, lv);
+ horizontal_bc(E, E->VB, top, 1, E->control.VBXtopval, SBX, 1, lv);
+ horizontal_bc(E, E->VB, top, 3, 0.0, SBZ, 0, lv);
+ horizontal_bc(E, E->VB, top, 2, E->control.VBYtopval, SBY, 1, lv);
}
}
@@ -76,24 +76,25 @@
for(lv = E->mesh.levmax; lv >= E->mesh.levmin; lv--)
{
- if(E->mesh.botvbc == 1)
+ top = E->mesh.NOZ[lv];
+ if(E->mesh.botvbc == 1) /* bottom velocity boundary condition */
{
- horizontal_bc(E, E->VB, 1, 1, E->control.VBXbotval, VBX, 1, lv);
- horizontal_bc(E, E->VB, 1, 3, 0.0, VBZ, 1, lv);
- horizontal_bc(E, E->VB, 1, 2, E->control.VBYbotval, VBY, 1, lv);
- horizontal_bc(E, E->VB, 1, 1, 0.0, SBX, 0, lv);
- horizontal_bc(E, E->VB, 1, 3, 0.0, SBZ, 0, lv);
- horizontal_bc(E, E->VB, 1, 2, 0.0, SBY, 0, lv);
+ horizontal_bc(E, E->VB, bottom, 1, E->control.VBXbotval, VBX, 1, lv);
+ horizontal_bc(E, E->VB, bottom, 3, 0.0, VBZ, 1, lv);
+ horizontal_bc(E, E->VB, bottom, 2, E->control.VBYbotval, VBY, 1, lv);
+ horizontal_bc(E, E->VB, bottom, 1, 0.0, SBX, 0, lv);
+ horizontal_bc(E, E->VB, bottom, 3, 0.0, SBZ, 0, lv);
+ horizontal_bc(E, E->VB, bottom, 2, 0.0, SBY, 0, lv);
}
- if(E->mesh.topvbc == 1)
+ if(E->mesh.topvbc == 1) /* top velocity boundary condition */
{
- E->control.VBXtopval = E->control.plate_vel;
- horizontal_bc(E, E->VB, E->mesh.NOZ[lv], 1, E->control.VBXtopval, VBX, 1, lv);
- horizontal_bc(E, E->VB, E->mesh.NOZ[lv], 3, 0.0, VBZ, 1, lv);
- horizontal_bc(E, E->VB, E->mesh.NOZ[lv], 2, E->control.VBYtopval, VBY, 1, lv);
- horizontal_bc(E, E->VB, E->mesh.NOZ[lv], 1, 0.0, SBX, 0, lv);
- horizontal_bc(E, E->VB, E->mesh.NOZ[lv], 3, 0.0, SBZ, 0, lv);
- horizontal_bc(E, E->VB, E->mesh.NOZ[lv], 2, 0.0, SBY, 0, lv);
+ E->control.VBXtopval = E->control.plate_vel; /* this should be more explicit */
+ horizontal_bc(E, E->VB, top, 1, E->control.VBXtopval, VBX, 1, lv);
+ horizontal_bc(E, E->VB, top, 3, 0.0, VBZ, 1, lv);
+ horizontal_bc(E, E->VB, top, 2, E->control.VBYtopval, VBY, 1, lv);
+ horizontal_bc(E, E->VB, top, 1, 0.0, SBX, 0, lv);
+ horizontal_bc(E, E->VB, top, 3, 0.0, SBZ, 0, lv);
+ horizontal_bc(E, E->VB, top, 2, 0.0, SBY, 0, lv);
}
}
@@ -535,28 +536,25 @@
void velocities_conform_bcs(struct All_variables *E, double *U)
{
- //int node, d;
- int node;
+ int node;
+
+ const unsigned int typex = VBX;
+ const unsigned int typez = VBZ;
+ const unsigned int typey = VBY;
+
+ const int nno = E->lmesh.nno;
- const unsigned int typex = VBX;
- const unsigned int typez = VBZ;
- const unsigned int typey = VBY;
-
- //const int dofs = E->mesh.dof;
- const int nno = E->lmesh.nno;
-
- for(node = 1; node <= nno; node++)
- {
- if(E->node[node] & typex)
- U[E->id[node].doff[1]] = E->VB[1][node];
- if(E->node[node] & typez)
- U[E->id[node].doff[3]] = E->VB[3][node];
- if(E->node[node] & typey)
- U[E->id[node].doff[2]] = E->VB[2][node];
-
- }
-
- return;
+ for(node = 1; node <= nno; node++){
+ if(E->node[node] & typex)
+ U[E->id[node].doff[1]] = E->VB[1][node];
+ if(E->node[node] & typey)
+ U[E->id[node].doff[2]] = E->VB[2][node];
+ if(E->node[node] & typez){
+ U[E->id[node].doff[3]] = E->VB[3][node];
+ }
+ }
+
+ return;
}
Modified: mc/3D/CitcomCU/trunk/src/Citcom.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Citcom.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Citcom.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -108,16 +108,18 @@
while(E.control.keep_going && (Emergency_stop == 0))
{
+ //if(E.parallel.me == 0)fprintf(stderr,"processing heating\n");
process_heating(&E);
E.monitor.solution_cycles++;
if(E.monitor.solution_cycles > E.control.print_convergence)
E.control.print_convergence = 1;
-
+ // if(E.parallel.me == 0)fprintf(stderr,"processing buoyancy\n");
/**/ report(&E, "Update buoyancy for further `timesteps'");
(E.next_buoyancy_field) (&E);
-
+ //if(E.parallel.me == 0)fprintf(stderr,"processing temp\n");
/**/ report(&E, "Process results of buoyancy update");
process_temp_field(&E, E.monitor.solution_cycles);
+ //if(E.parallel.me == 0)fprintf(stderr,"solving stokes\n");
general_stokes_solver(&E);
Modified: mc/3D/CitcomCU/trunk/src/Composition_adv.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Composition_adv.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Composition_adv.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -109,10 +109,8 @@
{
int i;
double temp1, temp2, temp3;
-
/* velocity VO at t and x=XMC */
velocity_markers(E, V, on_off);
-
/* predicted marker positions at t+dt */
if(E->control.CART3D)
{
@@ -135,13 +133,10 @@
E->XMCpred[3][i] = temp3;
}
}
-
transfer_markers_processors(E, on_off);
-
/* predicted compositional field at t+dt */
element_markers(E, on_off);
get_C_from_markers(E, C);
-
return;
}
@@ -151,7 +146,7 @@
{
//FILE *fp;
//char output_file[255];
- int i, proc, neighbor, no_transferred, no_received;
+ int i, proc, neighbor, no_transferred, no_received,asize;
static int been = 0;
static int markers;
@@ -162,23 +157,23 @@
if(been == 0)
{
markers = E->advection.markers / 10;
+ asize = (markers + 1) * E->mesh.nsd * 2;
for(neighbor = 1; neighbor <= E->parallel.no_neighbors; neighbor++)
{
E->parallel.traces_transfer_index[neighbor] = (int *)malloc((markers + 1) * sizeof(int));
- E->RVV[neighbor] = (float *)malloc((markers + 1) * E->mesh.nsd * 2 * sizeof(int));
- E->RXX[neighbor] = (double *)malloc((markers + 1) * E->mesh.nsd * 2 * sizeof(double));
+ E->RVV[neighbor] = (float *)malloc(asize * sizeof(int));
+ E->RXX[neighbor] = (double *)malloc(asize * sizeof(double));
E->RINS[neighbor] = (int *)malloc((markers + 1) * 2 * sizeof(int));
- E->PVV[neighbor] = (float *)malloc((markers + 1) * E->mesh.nsd * 2 * sizeof(int));
- E->PXX[neighbor] = (double *)malloc((markers + 1) * E->mesh.nsd * 2 * sizeof(double));
+ E->PVV[neighbor] = (float *)malloc(asize * sizeof(int));
+ E->PXX[neighbor] = (double *)malloc(asize * sizeof(double));
E->PINS[neighbor] = (int *)malloc((markers + 1) * 2 * sizeof(int));
}
E->traces_leave_index = (int *)malloc((markers + 1) * sizeof(int));
been++;
}
- for(neighbor = 1; neighbor <= E->parallel.no_neighbors; neighbor++)
+ for(neighbor = 0; neighbor <= E->parallel.no_neighbors; neighbor++)
E->parallel.traces_transfer_number[neighbor] = 0;
-
if(on_off == 1)
{ // use XMC
for(i = 1; i <= E->advection.markers; i++)
@@ -238,15 +233,13 @@
}
// prepare for transfer
-
+
no_transferred = 0;
for(neighbor = 1; neighbor <= E->parallel.no_neighbors; neighbor++)
{
no_transferred += E->parallel.traces_transfer_number[neighbor];
}
-
prepare_transfer_arrays(E);
-
exchange_number_rec_markers(E);
no_received = 0;
@@ -401,31 +394,38 @@
void prepare_transfer_arrays(struct All_variables *E)
{
- int j, i, neighbor, k1, k2, k3;
-
+ int j, part, neighbor, k1, k2, k3;
+ //if(E->parallel.me==0)fprintf(stderr,"ta 1 ok\n");
for(neighbor = 1; neighbor <= E->parallel.no_neighbors; neighbor++)
{
k1 = k2 = k3 = 0;
+ if((E->parallel.me==0) && (E->monitor.solution_cycles>199))
+ //fprintf(stderr,"ta %i %i %i - %i %i - %i %i %i \n",neighbor, E->parallel.no_neighbors,E->parallel.traces_transfer_number[neighbor],E->parallel.traces_transfer_number[neighbor]*6,E->parallel.traces_transfer_number[neighbor]*2,E->advection.markers / 10 ,(E->advection.markers / 10 + 1) * E->mesh.nsd * 2 ,(E->advection.markers / 10 + 1) *2);
for(j = 0; j < E->parallel.traces_transfer_number[neighbor]; j++)
{
- i = E->parallel.traces_transfer_index[neighbor][j];
- E->PVV[neighbor][k1++] = E->VO[1][i];
- E->PVV[neighbor][k1++] = E->VO[2][i];
- E->PVV[neighbor][k1++] = E->VO[3][i];
- E->PVV[neighbor][k1++] = E->Vpred[1][i];
- E->PVV[neighbor][k1++] = E->Vpred[2][i];
- E->PVV[neighbor][k1++] = E->Vpred[3][i];
- E->PXX[neighbor][k2++] = E->XMC[1][i];
- E->PXX[neighbor][k2++] = E->XMC[2][i];
- E->PXX[neighbor][k2++] = E->XMC[3][i];
- E->PXX[neighbor][k2++] = E->XMCpred[1][i];
- E->PXX[neighbor][k2++] = E->XMCpred[2][i];
- E->PXX[neighbor][k2++] = E->XMCpred[3][i];
- E->PINS[neighbor][k3++] = E->C12[i];
- E->PINS[neighbor][k3++] = E->CElement[i];
+ part = E->parallel.traces_transfer_index[neighbor][j];
+ if((part > E->advection.markers)||(part<1)){fprintf(stderr,"pta: out of bounds %i %i\n",part,E->advection.markers);}
+ if((E->parallel.me==0) && (E->monitor.solution_cycles>199))fprintf(stderr,"%i %i %i\n",neighbor,j,part);
+ E->PVV[neighbor][k1++] = E->VO[1][part];
+ E->PVV[neighbor][k1++] = E->VO[2][part];
+ E->PVV[neighbor][k1++] = E->VO[3][part];
+ E->PVV[neighbor][k1++] = E->Vpred[1][part];
+ E->PVV[neighbor][k1++] = E->Vpred[2][part];
+ E->PVV[neighbor][k1++] = E->Vpred[3][part];
+
+ E->PXX[neighbor][k2++] = E->XMC[1][part];
+ E->PXX[neighbor][k2++] = E->XMC[2][part];
+ E->PXX[neighbor][k2++] = E->XMC[3][part];
+ E->PXX[neighbor][k2++] = E->XMCpred[1][part];
+ E->PXX[neighbor][k2++] = E->XMCpred[2][part];
+ E->PXX[neighbor][k2++] = E->XMCpred[3][part];
+
+ E->PINS[neighbor][k3++] = E->C12[part];
+ E->PINS[neighbor][k3++] = E->CElement[part];
}
+ //if((E->parallel.me==0) && (E->monitor.solution_cycles>199))fprintf(stderr,"ta inner loop ok\n");
}
-
+ //if(E->parallel.me==0)fprintf(stderr,"ta 2 ok\n");
return;
}
Modified: mc/3D/CitcomCU/trunk/src/Construct_arrays.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Construct_arrays.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Construct_arrays.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -380,6 +380,7 @@
const double zero = 0.0;
+
for(level = E->mesh.levmax; level >= E->mesh.levmin; level--)
{
neq = E->lmesh.NEQ[level];
@@ -528,7 +529,7 @@
}
}
- }
+ } /* level loop */
return;
}
@@ -666,72 +667,66 @@
void construct_elt_ks(struct All_variables *E)
{
- //int e, el, lev, j, k, ii;
- int el, lev, j, k, ii;
-
- const int dims = E->mesh.nsd;
- const int n = loc_mat_size[E->mesh.nsd];
-
- if(E->control.verbose && E->parallel.me == 0)
- fprintf(stderr, "storing elt k matrices\n");
- if(E->parallel.me == 0)
- fprintf(stderr, "storing elt k matrices\n");
-
- for(lev = E->mesh.levmin; lev <= E->mesh.levmax; lev++)
+ //int e, el, lev, j, k, ii;
+ int el, lev, j, k, ii;
+ const int dims = E->mesh.nsd;
+ const int n = loc_mat_size[E->mesh.nsd];
+
+ if(E->control.verbose && E->parallel.me == 0)
+ fprintf(stderr, "storing elt k matrices\n");
+
+ for(lev = E->mesh.levmin; lev <= E->mesh.levmax; lev++)
+ {
+
+ E->parallel.idb = 1;
+
+ for(el = 1; el <= E->lmesh.NEL[lev]; el++)
{
-
- E->parallel.idb = 1;
-
- for(el = 1; el <= E->lmesh.NEL[lev]; el++)
- {
-
- get_elt_k(E, el, E->elt_k[lev][el].k, lev, 0); /* not for penalty */
-
- if(E->control.augmented_Lagr)
- get_aug_k(E, el, E->elt_k[lev][el].k, lev);
-
- build_diagonal_of_K(E, el, E->elt_k[lev][el].k, lev);
-
-
- }
-
- exchange_id_d20(E, E->BI[lev], lev); /*correct BI */
-
- for(j = 0; j < E->lmesh.NEQ[lev]; j++)
- {
- if(E->BI[lev][j] == 0.0)
- fprintf(stderr, "me= %d level %d, equation %d/%d has zero diagonal term\n", E->parallel.me, lev, j, E->mesh.NEQ[lev]);
- assert(E->BI[lev][j] != 0 /* diagonal of matrix = 0, not acceptable */ );
- E->BI[lev][j] = (float)1.0 / E->BI[lev][j];
- }
+
+ get_elt_k(E, el, E->elt_k[lev][el].k, lev, 0); /* not for penalty */
+
+ if(E->control.augmented_Lagr)
+ get_aug_k(E, el, E->elt_k[lev][el].k, lev);
+
+ build_diagonal_of_K(E, el, E->elt_k[lev][el].k, lev);
+
+
}
-
- if(E->control.verbose)
- for(lev = E->mesh.levmin; lev <= E->mesh.levmax; lev++)
- for(el = 1; el <= E->lmesh.NEL[lev]; el++)
- for(j = 1; j <= enodes[E->mesh.nsd]; j++)
- for(k = 1; k <= enodes[E->mesh.nsd]; k++)
- {
- ii = (j * n + k) * dims - (dims * n + dims);
- /* fprintf(E->fp,"stiff_for_e %d %d %d %g %g %g %g \n",el,j,k,E->elt_k[lev][el].k[ii],E->elt_k[lev][el].k[ii+1],E->elt_k[lev][el].k[ii+n],E->elt_k[lev][el].k[ii+n+1]); */
- }
-
- return;
+
+ exchange_id_d20(E, E->BI[lev], lev); /*correct BI */
+
+ for(j = 0; j < E->lmesh.NEQ[lev]; j++)
+ {
+ if(E->BI[lev][j] == 0.0)
+ fprintf(stderr, "me= %d level %d, equation %d/%d has zero diagonal term\n", E->parallel.me, lev, j, E->mesh.NEQ[lev]);
+ assert(E->BI[lev][j] != 0 /* diagonal of matrix = 0, not acceptable */ );
+ E->BI[lev][j] = (float)1.0 / E->BI[lev][j];
+ }
+ }
+
+
+ if(E->control.verbose)
+ for(lev = E->mesh.levmin; lev <= E->mesh.levmax; lev++){
+ for(el = 1; el <= E->lmesh.NEL[lev]; el++)
+ for(j = 1; j <= enodes[E->mesh.nsd]; j++)
+ for(k = 1; k <= enodes[E->mesh.nsd]; k++)
+ {
+ ii = (j * n + k) * dims - (dims * n + dims);
+ /* fprintf(E->fp,"stiff_for_e %d %d %d %g %g %g %g \n",el,j,k,E->elt_k[lev][el].k[ii],E->elt_k[lev][el].k[ii+1],E->elt_k[lev][el].k[ii+n],E->elt_k[lev][el].k[ii+n+1]); */
+ }
+ }
+
+ return;
}
void construct_elt_gs(struct All_variables *E)
{
- //int el, lev, a;
int el, lev;
- //const int dims = E->mesh.nsd;
- //const int dofs = E->mesh.dof;
- //const int ends = enodes[dims];
-
if(E->control.verbose && E->parallel.me == 0)
- fprintf(stderr, "storing elt g matrices\n");
+ fprintf(stderr, "storing elt g matrices\n");
for(lev = E->mesh.levmin; lev <= E->mesh.levmax; lev++)
for(el = 1; el <= E->lmesh.NEL[lev]; el++)
@@ -851,7 +846,6 @@
E->elt_k[i] = (struct EK *)malloc((E->lmesh.NEL[i] + 1) * sizeof(struct EK));
}
}
-
if(been_here0 == 0 || (E->viscosity.update_allowed && E->monitor.solution_cycles % E->control.KERNEL == 0))
{
/* do the following for the 1st time or update_allowed is true */
Modified: mc/3D/CitcomCU/trunk/src/Drive_solvers.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Drive_solvers.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Drive_solvers.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -80,11 +80,12 @@
damp = 0;
}
} /* end iterate */
- oldU = (double *)malloc(neq * sizeof(double));
- for(i = 0; i < neq; i++)
- oldU[i] = 0.0;
- visits++;
+ oldU = (double *)malloc(neq * sizeof(double));
+ for(i = 0; i < neq; i++)
+ oldU[i] = 0.0;
+ visits++;
}
+ //if(E->parallel.me==0)fprintf(stderr,"Stoked prep done\n");
dUdot_mag = 0.0;
@@ -99,7 +100,9 @@
time = CPU_time0();
velocities_conform_bcs(E, E->U);
+ //if(E->parallel.me==0)fprintf(stderr,"assembling forces\n");
assemble_forces(E, 0);
+
/*
if(E->parallel.me==0)
@@ -114,8 +117,9 @@
do
{
if(E->viscosity.update_allowed)
- get_system_viscosity(E, 1, E->EVI[E->mesh.levmax], E->VI[E->mesh.levmax]);
+ get_system_viscosity(E, 1, E->EVI[E->mesh.levmax], E->VI[E->mesh.levmax]);
construct_stiffness_B_matrix(E);
+ //if(E->parallel.me==0)fprintf(stderr,"calling solver\n");
solve_constrained_flow_iterative(E);
E->monitor.visc_iter_count++;
@@ -145,28 +149,32 @@
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);
}
- E->monitor.visc_iter_count++;
- } /* end for SDEPV / BDEPV */
- } while((E->monitor.visc_iter_count < 50) && (dUdot_mag > E->viscosity.sdepv_misfit) && iterate);
+ } /* end for stress type iterations */
+ } while(iterate &&
+ (dUdot_mag > E->viscosity.sdepv_misfit) &&
+ (E->monitor.visc_iter_count < 50) );
+
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
+ if(E->viscosity.anivisc_start_from_iso) /* first step will be
solved isotropically at
first */
- return TRUE;
+ return 1;
else
- return (E->viscosity.SDEPV || E->viscosity.BDEPV)?(TRUE):(FALSE);
+ return (E->viscosity.SDEPV || E->viscosity.BDEPV)?(1):(0);
+ }else{
+#endif
+ /* regular operation */
+ return ((E->viscosity.SDEPV || E->viscosity.BDEPV)?(1):(0));
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
}
-#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-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Element_calculations.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -116,8 +116,9 @@
fprintf(E->fp,"bb %d %g\n",a, E->F[a]);
*/
+ //if(E->parallel.me == 0)fprintf(stderr,"af exchange\n");
exchange_id_d20(E, E->F, lev);
-
+ //if(E->parallel.me == 0)fprintf(stderr,"af strip BC\n");
strip_bcs_from_residual(E, E->F, lev);
return;
@@ -152,18 +153,14 @@
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;
if(E->control.Rsphere) /* need rtf for spherical */
- get_rtf(E, el, 0, rtf, lev);
+ get_rtf(E, el, 0, rtf, lev); /* vpts */
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];
@@ -178,12 +175,12 @@
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);
+ get_ba(&(E->N),&(E->GNX[lev][el]),&Cc, &Ccx,rtf, E->mesh.nsd,vpts, ends,1,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);
+ get_ba(&(E->N),&(E->GNX[lev][el]),&Cc, &Ccx,rtf, E->mesh.nsd, vpts, ends,0,ba);
}
#endif
; /* only anisotropic cartesian needs ba factors */
@@ -198,11 +195,15 @@
if(E->viscosity.allow_anisotropic_viscosity){
for(i = 1; i <= dims; i++)
for(j = 1; j <= dims; j++)
- for(k = 1; k <= VPOINTS3D; k++){
+ for(k = 1; k <= vpts; 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] );
+ 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]+
@@ -218,7 +219,7 @@
*/
if(E->control.CART3D){
/* cartesian isotropic does not use ba[] */
- for(k = 1; k <= VPOINTS3D; k++){
+ for(k = 1; k <= vpts; 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)];
@@ -231,7 +232,7 @@
}
temp = 0.0;
- for(k = 1; k <= VPOINTS3D; k++){
+ for(k = 1; k <= vpts; 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)]);
@@ -245,7 +246,7 @@
}else if(E->control.Rsphere){
for(i = 1; i <= dims; i++)
for(j = 1; j <= dims; j++)
- for(k = 1; k <= VPOINTS3D; k++){
+ for(k = 1; k <= vpts; 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]) +
@@ -289,57 +290,64 @@
}
/*
-compute the displacement gradient matrix B
+compute the displacement gradient matrix B at the velocity points
+(for element K computations)
*/
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])
+ int dims,int vpts, int ends, int spherical, double ba[9][4][9][7])
{
- int i,k,a,vpts,ends;
+ int i,k,a;
double ra[9], si[9], ct[9];
const double one = 1.0;
const double two = 2.0;
- float gnx0, gnx1, gnx2;
+ double gnx0, gnx1, gnx2;
double shp, cc1, cc2, cc3;
- vpts = vpoints[dims];
- ends = enodes[dims];
-
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];
+ ra[k] = rtf[3][k]; /* 1/r */
+ si[k] = one / sin(rtf[1][k]); /* 1/sin(t) */
+ ct[k] = cos(rtf[1][k]) * si[k]; /* 1/tan(t) */
}
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];
-
+ for(k = 1; k <= vpts; k++){
+ gnx0 = GNx->vpt[GNVXINDEX(0, a, k)]; /* d_t */
+ gnx1 = GNx->vpt[GNVXINDEX(1, a, k)]; /* d_p */
+ gnx2 = GNx->vpt[GNVXINDEX(2, a, k)]; /* d_r */
+ shp = N->vpt[GNVINDEX(a, k)];
+ for(i = 1; i <= dims; i++){
+ cc1 = cc->vpt[BVINDEX(1, i, a, k)]; /* t */
+ cc2 = cc->vpt[BVINDEX(2, i, a, k)]; /* p */
+ cc3 = cc->vpt[BVINDEX(3, i, a, k)]; /* r */
+ /* tt = (d_t ut + ur)/r -- 11 */
+ ba[a][i][k][1] = ((gnx0 * cc1 + shp * ccx->vpt[BVXINDEX(1, i, 1, a, k)]) +
+ shp * cc3) * ra[k];
+ /* pp = (ut/tan(t) + ur + (d_p up)/sin(t))/r -- 22 */
+ 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];
+ /* rr = d_r ur -- 33 */
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];
+ /* tp = (d_t up + up/tan(t) + d_p ut/sin(t) ) / r -- 21 + 12*/
+ 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];
+ /* tr = d_r ut + (d_t ur - ut)/r -- 13 + 31 */
+ ba[a][i][k][5] = gnx2 * cc1 +
+ (gnx0 * cc3 +
+ shp * (ccx->vpt[BVXINDEX(3, i, 1, a, k)] - cc1)) * ra[k];
+ /* pr = d_r up - up/r - d_p ur/sin(t)/r -- 23 + 32 */
+ ba[a][i][k][6] = gnx2 * cc2 +
+ (( gnx1 * cc3 + shp * ccx->vpt[BVXINDEX(3,i,2,a,k)] ) * si[k] - shp * cc2 ) * ra[k];
}
+ }
}else{ /* cartesian */
for(a = 1; a <= ends; a++)
- for(k = 1; k <= VPOINTS3D; k++){
+ for(k = 1; k <= vpts; k++){
gnx0 = GNx->vpt[GNVXINDEX(0, a, k)];
gnx1 = GNx->vpt[GNVXINDEX(1, a, k)];
gnx2 = GNx->vpt[GNVXINDEX(2, a, k)];
@@ -369,16 +377,235 @@
ba[a][3][k][6] = gnx1;
}
- } /* end Caretsian */
+ } /* end Cartesian */
}
+/*
- /* =============================================
- * General calling function for del_squared:
- * according to whether it should be element by
- * element or node by node.
- * ============================================= */
+ get B at the pressure integration points for strain-rate
+ computations
+ */
+void get_ba_p(struct Shape_function *N,
+ struct Shape_function_dx *GNx,
+ struct CC *cc, struct CCX *ccx, double rtf[4][9],
+ int dims,int ppts, int ends, int spherical,
+ double ba[9][4][9][7])
+{
+
+ int i,k,a;
+ double ra[9], si[9], ct[9];
+ const double one = 1.0;
+ const double two = 2.0;
+ double gnx0, gnx1, gnx2;
+ double shp, cc1, cc2, cc3;
+ if(spherical){ /* spherical coordinates */
+ for(k = 1; k <= ppts; k++){
+ ra[k] = rtf[3][k];
+ si[k] = one / sin(rtf[1][k]);
+ ct[k] = cos(rtf[1][k]) * si[k]; /* cot(t) */
+ }
+ for(a = 1; a <= ends; a++)
+ for(k = 1; k <= ppts; k++){
+ gnx0 = GNx->ppt[GNPXINDEX(0, a, k)];
+ gnx1 = GNx->ppt[GNPXINDEX(1, a, k)];
+ gnx2 = GNx->ppt[GNPXINDEX(2, a, k)];
+ shp = N->ppt[GNPINDEX(a, k)];
+ for(i = 1; i <= dims; i++){
+ cc1 = cc->ppt[BPINDEX(1, i, a, k)];
+ cc2 = cc->ppt[BPINDEX(2, i, a, k)];
+ cc3 = cc->ppt[BPINDEX(3, i, a, k)];
+
+ /* d_tt */
+ ba[a][i][k][1] = (gnx0 * cc1 + shp * ccx->ppt[BPXINDEX(1, i, 1, a, k)] + shp * cc3) * ra[k];
+ /* d_pp */
+ ba[a][i][k][2] = (shp * cc1 * ct[k] + shp * cc3 + (gnx1 * cc2 + shp * ccx->ppt[BPXINDEX(2, i, 2, a, k)]) * si[k]) * ra[k];
+ /* d_rr */
+ ba[a][i][k][3] = gnx2 * cc3;
+ /* d_tp */
+ ba[a][i][k][4] = (gnx0 * cc2 + shp * ccx->ppt[BPXINDEX(2, i, 1, a, k)] - shp * cc2 * ct[k] + (gnx1 * cc1 + shp * ccx->ppt[BPXINDEX(1, i, 2, a, k)]) * si[k]) * ra[k];
+ /* d_tr */
+ ba[a][i][k][5] = gnx2 * cc1 + (gnx0 * cc3 + shp * (ccx->ppt[BPXINDEX(3, i, 1, a, k)] - cc1)) * ra[k];
+ /* d_pr */
+ ba[a][i][k][6] = gnx2 * cc2 + ((gnx1 * cc3 + shp * ccx->ppt[BPXINDEX(3, i, 2, a, k)]) * si[k] - shp * cc2 ) * ra[k];
+ }
+ }
+ }else{
+ /* cartesian coordinates */
+ for(a = 1; a <= ends; a++)
+ for(k = 1; k <= ppts; k++){
+ gnx0 = GNx->ppt[GNPXINDEX(0, a, k)];
+ gnx1 = GNx->ppt[GNPXINDEX(1, a, k)];
+ gnx2 = GNx->ppt[GNPXINDEX(2, a, k)];
+ /* xx */
+ ba[a][1][k][1] = gnx0;
+ ba[a][2][k][1] = 0.;
+ ba[a][3][k][1] = 0.;
+ /* yy */
+ ba[a][1][k][2] = 0.;
+ ba[a][2][k][2] = gnx1;
+ ba[a][3][k][2] = 0.;
+ /* zz */
+ ba[a][1][k][3] = 0.;
+ ba[a][2][k][3] = 0.;
+ ba[a][3][k][3] = gnx2;
+ /* xy */
+ ba[a][1][k][4] = gnx1;
+ ba[a][2][k][4] = gnx0;
+ ba[a][3][k][4] = 0.;
+ /* xz */
+ ba[a][1][k][5] = gnx2;
+ ba[a][2][k][5] = 0.;
+ ba[a][3][k][5] = gnx0;
+ /* yz */
+ ba[a][1][k][6] = 0.;
+ ba[a][2][k][6] = gnx2;
+ ba[a][3][k][6] = gnx1;
+
+ }
+ } /* end Cartesian */
+}
+
+/*
+
+ get velocity gradient matrix at element, and also compute the
+ average velocity in this element
+
+
+*/
+
+void get_vgm_p(double VV[4][9],struct Shape_function *N,
+ struct Shape_function_dx *GNx,
+ struct CC *cc, struct CCX *ccx, double rtf[4][9],
+ int dims,int ppts, int ends, int spherical,
+ double l[3][3], double v[3])
+{
+
+ int i,k,j,a;
+ double ra[9], si[9], ct[9];
+ const double one = 1.0;
+ const double two = 2.0;
+ double vgm[3][3];
+ double shp, cc1, cc2, cc3,d_t,d_r,d_p,up,ur,ut;
+ /* init L matrix */
+ for(i=0;i < 3;i++){
+ v[i] = 0.0;
+ for(j=0;j < 3; j++)
+ l[i][j] = 0.0;
+ }
+ /* mean velocity at pressure integration point */
+ for(a=1;a <= ends;a++){
+ v[0] += N->ppt[GNPINDEX(a, 1)] * VV[1][a];
+ v[1] += N->ppt[GNPINDEX(a, 1)] * VV[2][a];
+ v[2] += N->ppt[GNPINDEX(a, 1)] * VV[3][a];
+ }
+ if(spherical){
+ for(k = 1; k <= ppts; k++){
+ ra[k] = rtf[3][k]; /* 1/r */
+ si[k] = one / sin(rtf[1][k]); /* 1/sin(t) */
+ ct[k] = cos(rtf[1][k]) * si[k]; /* cot(t) */
+ }
+ for(a = 1; a <= ends; a++){
+ for(k = 1; k <= ppts; k++){
+ d_t = GNx->ppt[GNPXINDEX(0, a, k)]; /* d_t */
+ d_p = GNx->ppt[GNPXINDEX(1, a, k)]; /* d_p */
+ d_r = GNx->ppt[GNPXINDEX(2, a, k)]; /* d_r */
+ shp = N->ppt[GNPINDEX(a, k)];
+ for(i = 1; i <= dims; i++){
+ ut = cc->ppt[BPINDEX(1, i, a, k)]; /* ut */
+ up = cc->ppt[BPINDEX(2, i, a, k)]; /* up */
+ ur = cc->ppt[BPINDEX(3, i, a, k)]; /* ur */
+
+ /* velocity gradient matrix is transpose of grad v, using Citcom sort t, p, r
+
+ | d_t(vt) d_p(vt) d_r(vt) |
+ | d_t(vp) d_p(vp) d_r(vp) |
+ | d_t(vr) d_p(vr) d_r(vr) |
+
+ */
+
+ /* d_t vt = 1/r (d_t vt + vr) */
+ vgm[0][0] = ((d_t * ut + shp * ccx->ppt[BPXINDEX(1, i, 1, a, k)]) +
+ shp * ur) * ra[k];
+ /* d_p vt = 1/r (1/sin(t) d_p vt -vp/tan(t)) */
+ vgm[0][1] = ((d_p * ut + shp * ccx->ppt[BPXINDEX(1, i, 2, a, k)]) * si[k] -
+ shp * up * ct[k]) * ra[k];
+ /* d_r vt = d_r v_t */
+ vgm[0][2] = d_r * ut;
+ /* d_t vp = 1/r d_t v_p*/
+ vgm[1][0] = (d_t * up + shp * ccx->ppt[BPXINDEX(2, i, 1, a, k)]) * ra[k];
+ /* d_p vp = 1/r((d_p vp)/sin(t) + vt/tan(t) + vr) */
+ vgm[1][1] = ((d_p * up + shp * ccx->ppt[BPXINDEX(2, i, 2, a, k)]) * si[k] +
+ shp * ut * ct[k] + shp * ur) * ra[k];
+ /* d_r vp = d_r v_p */
+ vgm[1][2] = d_r * up;
+ /* d_t vr = 1/r(d_t vr - vt) */
+ vgm[2][0] = ((d_t * ur + shp * ccx->ppt[BPXINDEX(3, i, 1, a, k)]) -
+ shp * ut) * ra[k];
+ /* d_p vr = 1/r(1/sin(t) d_p vr - vp) */
+ vgm[2][1] = (( d_p * ur + shp * ccx->ppt[BPXINDEX(3,i, 2,a,k)] ) * si[k] -
+ shp * up ) * ra[k];
+ /* d_r vr = d_r vr */
+ vgm[2][2] = d_r * ur;
+
+
+ l[0][0] += vgm[0][0] * VV[i][a];
+ l[0][1] += vgm[0][1] * VV[i][a];
+ l[0][2] += vgm[0][2] * VV[1][a];
+
+ l[1][0] += vgm[1][0] * VV[i][a];
+ l[1][1] += vgm[1][1] * VV[i][a];
+ l[1][2] += vgm[1][2] * VV[i][a];
+
+ l[2][0] += vgm[2][0] * VV[i][a];
+ l[2][1] += vgm[2][1] * VV[i][a];
+ l[2][2] += vgm[2][2] * VV[i][a];
+
+ }
+ }
+ }
+ }else{
+ /* cartesian */
+ for(k = 1; k <= ppts; k++){
+ for(a = 1; a <= ends; a++){
+ /* velocity gradient matrix is transpose of grad v
+
+ | d_x(vx) d_y(vx) d_z(vx) |
+ | d_x(vy) d_y(vy) d_z(vy) |
+ | d_x(vz) d_y(vz) d_z(vz) |
+ */
+ l[0][0] += GNx->ppt[GNPXINDEX(0, a, k)] * VV[1][a]; /* other contributions are zero */
+ l[0][1] += GNx->ppt[GNPXINDEX(1, a, k)] * VV[1][a];
+ l[0][2] += GNx->ppt[GNPXINDEX(2, a, k)] * VV[1][a];
+
+ l[1][0] += GNx->ppt[GNPXINDEX(0, a, k)] * VV[2][a];
+ l[1][1] += GNx->ppt[GNPXINDEX(1, a, k)] * VV[2][a];
+ l[1][2] += GNx->ppt[GNPXINDEX(2, a, k)] * VV[2][a];
+
+ l[2][0] += GNx->ppt[GNPXINDEX(0, a, k)] * VV[3][a];
+ l[2][1] += GNx->ppt[GNPXINDEX(1, a, k)] * VV[3][a];
+ l[2][2] += GNx->ppt[GNPXINDEX(2, a, k)] * VV[3][a];
+
+ }
+ }
+ }
+ if(ppts != 1){
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++)
+ l[i][j] /= (float)ppts;
+ }
+
+}
+
+
+
+
+/* =============================================
+ * General calling function for del_squared:
+ * according to whether it should be element by
+ * element or node by node.
+ * ============================================= */
+
void assemble_del2_u(struct All_variables *E, double *u, double *Au, int level, int strip_bcs)
{
if(E->control.NMULTIGRID || E->control.NASSEMBLE)
@@ -694,7 +921,7 @@
{
p = (a - 1) * dims;
j = E->LMD[level][e].node[a].doff[1];
- gradP[p] += E->BI[level][j] * E->elt_del[level][e].g[p][0];
+ gradP[p] += E->BI[level][j] * E->elt_del[level][e].g[p][0];
j = E->LMD[level][e].node[a].doff[2];
gradP[p + 1] += E->BI[level][j] * E->elt_del[level][e].g[p + 1][0];
@@ -732,70 +959,105 @@
void get_elt_g(struct All_variables *E, int el, higher_precision elt_del[24][1], int lev)
{
- //double dGNdash[3];
- //double recip_radius, temp;
- double temp;
- //int p, a, nint, es, d, i, j, k;
- int p, a, i;
- double ra, ct, si, x[4], rtf[4][9];
- //int lmsize;
+ double temp;
+ int p, a, i;
+ double ra, ct,si, x[4], rtf[4][9];
+ static struct CC Cc;
+ static struct CCX Ccx;
+ const int dims = E->mesh.nsd;
+ const int ends = enodes[dims];
+ const int ppts = ppoints[dims];
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ const int modify_g = 1;
+ //const int modify_g = 0;
+ int j,k,off;
+ double Dtmp[6][6],Duse[6][6],rtf2[4][9],weight;
+ double ba[9][4][9][7];
+ const int vpts = vpoints[dims];
+#endif
- //struct Shape_function GN;
- //struct Shape_function_dA dOmega;
- //struct Shape_function_dx GNx;
- static struct CC Cc;
- static struct CCX Ccx;
- const int dims = E->mesh.nsd;
- //const int dofs = E->mesh.dof;
- const int ends = enodes[dims];
- //const int vpts = vpoints[dims];
- //const int sphere_key = 1;
-
- /* Special case, 4/8 node bilinear cartesian square/cube element -> 1 pressure point */
-
- temp = p_point[1].weight[dims - 1] * E->GDA[lev][el].ppt[1];
-
- /* unroll etc some more */
-
- if(E->control.Rsphere)
- {
- if((el - 1) % E->lmesh.ELZ[lev] == 0)
- construct_c3x3matrix_el(E, el, &Cc, &Ccx, lev, 1);
-
- get_rtf(E, el, 2, rtf, lev);
-
- ra = rtf[3][1];
- si = 1.0 / sin(rtf[1][1]);
- ct = cos(rtf[1][1]) * si;
-
- for(a = 1; a <= ends; a++)
- {
- for(i = 1; i <= dims; i++)
- x[i] = E->GNX[lev][el].ppt[GNPXINDEX(2, a, 1)] * Cc.ppt[BPINDEX(3, i, a, 1)]
- + 2.0 * ra * E->N.ppt[GNPINDEX(a, 1)] * Cc.ppt[BPINDEX(3, i, a, 1)] + ra * (E->GNX[lev][el].ppt[GNPXINDEX(0, a, 1)] * Cc.ppt[BPINDEX(1, i, a, 1)] + E->N.ppt[GNPINDEX(a, 1)] * Ccx.ppt[BPXINDEX(1, i, 1, a, 1)] + ct * E->N.ppt[GNPINDEX(a, 1)] * Cc.ppt[BPINDEX(1, i, a, 1)] + si * (E->GNX[lev][el].ppt[GNPXINDEX(1, a, 1)] * Cc.ppt[BPINDEX(2, i, a, 1)] + E->N.ppt[GNPINDEX(a, 1)] * Ccx.ppt[BPXINDEX(2, i, 2, a, 1)]));
-
- p = dims * (a - 1);
- elt_del[p][0] = -x[1] * temp;
- elt_del[p + 1][0] = -x[2] * temp;
- elt_del[p + 2][0] = -x[3] * temp;
-
- }
-
+ /* Special case, 4/8 node bilinear cartesian square/cube element ->
+ 1 pressure point */
+
+ temp = p_point[1].weight[dims - 1] * E->GDA[lev][el].ppt[1];
+
+ /* unroll etc some more */
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ if(E->viscosity.allow_anisotropic_viscosity){
+ if(E->control.Rsphere)
+ get_rtf(E, el, 0, rtf2, lev); /* vpts */
+ /* find avg constitutive matrix */
+ for(i=0;i<6;i++)
+ for(j=0;j<6;j++)
+ Duse[i][j]=0.0;
+ weight = 1./(2.*vpts);
+ for(i=1;i <= vpts;i++){
+ off = (el-1)*vpts+i;
+ get_constitutive(Dtmp,lev,off,rtf2[1][i],rtf2[2][i],(E->control.Rsphere),E);
+ for(j=0;j<6;j++)
+ for(k=0;k<6;k++)
+ Duse[j][k] += Dtmp[j][k]*weight;
+ }
+ if(E->control.Rsphere){
+ if((el - 1) % E->lmesh.ELZ[lev] == 0)
+ construct_c3x3matrix_el(E, el, &Cc, &Ccx, lev, 1);
+ get_rtf(E, el, 1, rtf, lev);
+ }
+ get_ba_p(&(E->N),&(E->GNX[lev][el]),&Cc, &Ccx,rtf,
+ E->mesh.nsd,ppts,ends,(E->control.Rsphere),ba);
+ /* assume single pressure point */
+ for(a = 1; a <= ends; a++){
+ for(i = 1; i <= dims; i++){
+ x[i] = 0.0;
+ for(k=0;k < 6;k++){
+ x[i] += Duse[0][k] * ba[a][i][1][k+1];
+ x[i] += Duse[1][k] * ba[a][i][1][k+1];
+ x[i] += Duse[2][k] * ba[a][i][1][k+1];
}
- else if(E->control.CART3D)
- {
+ }
+
+ p = dims * (a - 1);
+ elt_del[p][0] = -x[1] * temp; /* contributions of each velocity to pressure */
+ elt_del[p + 1][0] = -x[2] * temp;
+ elt_del[p + 2][0] = -x[3] * temp;
+ }
+ }else{
+#endif
+ if(E->control.Rsphere){
+ if((el - 1) % E->lmesh.ELZ[lev] == 0)
+ construct_c3x3matrix_el(E, el, &Cc, &Ccx, lev, 1);
+
+ get_rtf(E, el, 1, rtf, lev);
+
+ ra = rtf[3][1];
+ si = 1.0 / sin(rtf[1][1]);
+ ct = cos(rtf[1][1]) * si;
- for(a = 1; a <= ends; a++)
- {
- p = dims * (a - 1);
- elt_del[p][0] = -E->GNX[lev][el].ppt[GNPXINDEX(0, a, 1)] * temp;
- elt_del[p + 1][0] = -E->GNX[lev][el].ppt[GNPXINDEX(1, a, 1)] * temp;
- elt_del[p + 2][0] = -E->GNX[lev][el].ppt[GNPXINDEX(2, a, 1)] * temp;
- }
- }
-
- return;
+ for(a = 1; a <= ends; a++){
+ for(i = 1; i <= dims; i++)
+ x[i] = E->GNX[lev][el].ppt[GNPXINDEX(2, a, 1)] * Cc.ppt[BPINDEX(3, i, a, 1)]
+ + 2.0 * ra * E->N.ppt[GNPINDEX(a, 1)] * Cc.ppt[BPINDEX(3, i, a, 1)] + ra * (E->GNX[lev][el].ppt[GNPXINDEX(0, a, 1)] * Cc.ppt[BPINDEX(1, i, a, 1)] + E->N.ppt[GNPINDEX(a, 1)] * Ccx.ppt[BPXINDEX(1, i, 1, a, 1)] + ct * E->N.ppt[GNPINDEX(a, 1)] * Cc.ppt[BPINDEX(1, i, a, 1)] + si * (E->GNX[lev][el].ppt[GNPXINDEX(1, a, 1)] * Cc.ppt[BPINDEX(2, i, a, 1)] + E->N.ppt[GNPINDEX(a, 1)] * Ccx.ppt[BPXINDEX(2, i, 2, a, 1)]));
+
+ p = dims * (a - 1);
+ elt_del[p][0] = -x[1] * temp;
+ elt_del[p + 1][0] = -x[2] * temp;
+ elt_del[p + 2][0] = -x[3] * temp;
+
+ }
+ }else if(E->control.CART3D){
+ for(a = 1; a <= ends; a++){
+ p = dims * (a - 1);
+ elt_del[p][0] = -E->GNX[lev][el].ppt[GNPXINDEX(0, a, 1)] * temp;
+ elt_del[p + 1][0] = -E->GNX[lev][el].ppt[GNPXINDEX(1, a, 1)] * temp;
+ elt_del[p + 2][0] = -E->GNX[lev][el].ppt[GNPXINDEX(2, a, 1)] * temp;
+ }
+ }
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ }
+#endif
+
+ return;
}
@@ -854,30 +1116,17 @@
void get_elt_f(struct All_variables *E, int el, double elt_f[24], int penalty, int bcs)
{
- //int aid, i, p, a, b, d, j, k, q, es;
int i, p, a, b, j, k, q;
- //int node[5], back_front, got_elt_k, nodea, nodeb;
int got_elt_k, nodea, nodeb;
unsigned int type[4];
- //static int been_here = 0;
- //double force[9], force_at_gs[9], stress[9];
double force[9], force_at_gs[9];
- //double vector[4], magnitude;
- //double tmp, rtf[4][9];
- //double rtf[4][9];
double elt_k[24 * 24];
- //struct Shape_function GN;
- //struct Shape_function_dA dOmega;
- //struct Shape_function_dx GNx;
- //struct Shape_function1 GM;
- //struct Shape_function1_dA dGammax;
static struct CC Cc;
static struct CCX Ccx;
const int dims = E->mesh.nsd;
- //const int dofs = E->mesh.dof;
const int n = loc_mat_size[dims];
const int ends = enodes[dims];
const int vpts = vpoints[dims];
@@ -950,7 +1199,6 @@
else if(E->control.Rsphere)
{
- //get_rtf(E, el, 0, rtf, E->mesh.levmax);
if((el - 1) % E->lmesh.elz == 0)
construct_c3x3matrix_el(E, el, &Cc, &Ccx, E->mesh.levmax, 0);
Modified: mc/3D/CitcomCU/trunk/src/General_matrix_functions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/General_matrix_functions.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/General_matrix_functions.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -365,7 +365,8 @@
Iterative solver also using multigrid ........
=========================================================== */
-int solve_del2_u(struct All_variables *E, double *d0, double *F, double acc, int high_lev, int ic)
+int solve_del2_u(struct All_variables *E, double *d0, double *F,
+ double acc, int high_lev)
{
static int been_here = 0;
static int up_heavy, down_heavy, v_steps_high;
@@ -481,7 +482,8 @@
}
}
- if((count > 0) && (residual > r0 * 2.0) || (fabs(residual - prior_residual) < acc * 0.1 && (residual > acc * 10.0)))
+ if((count > 0) && (residual > r0 * 2.0) ||
+ (fabs(residual - prior_residual) < acc * 0.1 && (residual > acc * 10.0)))
convergent = 0;
else
{
Modified: mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -51,7 +51,7 @@
void convection_initial_temperature_ggrd(struct All_variables *E)
{
- int ll, mm, i, j, k, p, node, ii;
+ int ll, mm, i, j, k, p, node, ii,slice,hit;
double temp, temp1, temp2, temp3, base, radius, radius2;
unsigned int type;
@@ -86,7 +86,7 @@
/* 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;
+ bot_t = (E->mesh.bottbc) ? E->control.TBCbotval : 0.0;
top_t = (E->mesh.toptbc) ? E->control.TBCtopval : 1.0;
for(i=1;i<=E->lmesh.nno;i++) {
@@ -129,13 +129,23 @@
/*
initialize the GMT grid files
*/
- if(E->control.slab_slice){ /* only slice of x - depth */
-
- if(ggrd_grdtrack_init_general(FALSE,E->control.ggrd.temp.gfile,
- char_dummy,"",E->control.ggrd.temp.d,
- (E->parallel.me==0),
- FALSE))
- myerror("grdtrack x-z init error",E);
+ if(E->control.ggrd_slab_slice){ /* only a few slices of x - depth */
+ if(E->control.ggrd_slab_slice == 1){
+ if(ggrd_grdtrack_init_general(FALSE,E->control.ggrd.temp.gfile,
+ char_dummy,"",E->control.ggrd_ss_grd,
+ (E->parallel.me==0),
+ FALSE))
+ myerror("grdtrack x-z init error",E);
+ }else{ /* several slab slices */
+ for(slice=0;slice<E->control.ggrd_slab_slice;slice++){
+ sprintf(pfile,"%s.%i.grd",E->control.ggrd.temp.gfile,slice+1);
+ if(ggrd_grdtrack_init_general(FALSE,pfile,
+ char_dummy,"",(E->control.ggrd_ss_grd+slice),
+ (E->parallel.me==0),
+ FALSE))
+ myerror("grdtrack x-z init error",E);
+ }
+ }
}else{ /* 3-D */
if(ggrd_grdtrack_init_general(TRUE,E->control.ggrd.temp.gfile,
E->control.ggrd.temp.dfile,"",
@@ -175,38 +185,39 @@
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){
+ if(E->control.ggrd_slab_slice){
/*
slab slice input
*/
- if(E->control.Rsphere) {
- if(E->SX[1][node] <= E->control.slab_theta_bound)
- /* spherical interpolation */
- ggrd_grdtrack_interpolate_xy((double)E->SX[2][node] * ONEEIGHTYOVERPI,
- (double)E->SX[1][node],
- E->control.ggrd.temp.d,&tadd,
- FALSE);
- else{
- if(E->SX[3][node] == E->segment.zzlayer[E->segment.zlayers-1])
- tadd = E->control.TBCtopval;
- else
- tadd = 1.0;
+ for(slice=hit=0;(!hit) && (slice < E->control.ggrd_slab_slice);slice++){
+ if(E->control.Rsphere) {
+ if(in_slab_slice(E->SX[1][node],slice,E)){
+ /* spherical interpolation */
+ ggrd_grdtrack_interpolate_xy((double)E->SX[2][node] * ONEEIGHTYOVERPI,
+ (double)E->SX[1][node],
+ (E->control.ggrd_ss_grd+slice),&tadd,
+ FALSE);
+ hit=1;
+ }
+ }else{ /* cartesian interpolation */
+ if(in_slab_slice(E->X[2][node],slice,E)){
+ ggrd_grdtrack_interpolate_xy((double)E->X[1][node],
+ (double)E->X[3][node],
+ (E->control.ggrd_ss_grd+slice),&tadd,
+ FALSE);
+ hit = 1;
+ }
}
- }else{ /* cartesian interpolation */
- if(E->X[2][node] <= E->control.slab_theta_bound)
- ggrd_grdtrack_interpolate_xy((double)E->X[1][node],
- (double)E->X[3][node],
- E->control.ggrd.temp.d,&tadd,
- FALSE);
- else{
- if(E->X[3][node] == E->segment.zzlayer[E->segment.zlayers-1])
- tadd = E->control.TBCtopval;
- else
- tadd = 1.0;
- }
}
+ if(!hit)
+ if(((E->control.Rsphere) && (E->SX[3][node] == E->segment.zzlayer[E->segment.zlayers-1]))||
+ (E->X[3][node] == E->segment.zzlayer[E->segment.zlayers-1]))
+ tadd = E->control.TBCtopval;
+ else
+ tadd = 1.0;
+ /* end slice part */
}else{
/*
@@ -335,7 +346,6 @@
y1 = E->SX[2][node];
/* linear */
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 */
@@ -346,7 +356,8 @@
ll = E->convection.perturb_ll[p];
con = E->convection.perturb_mag[p];
- E->T[node] += E->convection.perturb_mag[p] * modified_plgndr_a(ll, mm, x1) * cos(mm * y1) *
+ E->T[node] += E->convection.perturb_mag[p] *
+ modified_plgndr_a(ll, mm, x1) * cos(mm * y1) *
sin(M_PI * (z1 - E->sphere.ri) / (E->sphere.ro - E->sphere.ri));
}
}
@@ -371,14 +382,23 @@
mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), mpi_tag,
MPI_COMM_WORLD, &mpi_stat);
}
- if(E->control.slab_slice){
- /* x - z slice */
- if(ggrd_grdtrack_init_general(FALSE,E->control.ggrd.comp.gfile,
- char_dummy,"",
- E->control.ggrd.comp.d,
- /* (E->parallel.me==0)*/
- FALSE,FALSE))
- myerror("grdtrack init error",E);
+ if(E->control.ggrd_slab_slice){
+ if(E->control.ggrd_slab_slice == 1){
+ if(ggrd_grdtrack_init_general(FALSE,E->control.ggrd.comp.gfile,
+ char_dummy,"",E->control.ggrd_ss_grd,
+ (E->parallel.me==0),
+ FALSE))
+ myerror("grdtrack x-z init error",E);
+ }else{ /* several slab slices */
+ for(slice=0;slice<E->control.ggrd_slab_slice;slice++){
+ sprintf(pfile,"%s.%i.grd",E->control.ggrd.comp.gfile,slice+1);
+ if(ggrd_grdtrack_init_general(FALSE,pfile,
+ char_dummy,"",(E->control.ggrd_ss_grd+slice),
+ (E->parallel.me==0),
+ FALSE))
+ myerror("grdtrack x-z init error",E);
+ }
+ }
}else{
/* 3-D */
if(ggrd_grdtrack_init_general(TRUE,E->control.ggrd.comp.gfile,
@@ -401,27 +421,32 @@
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){
+ if(E->control.ggrd_slab_slice){
/* slab */
- if(E->control.Rsphere) {
- if(E->SX[1][node] <= E->control.slab_theta_bound)
- /* spherical interpolation */
- ggrd_grdtrack_interpolate_xy((double)E->SX[2][node] * ONEEIGHTYOVERPI,
- (double)E->SX[1][node],
- E->control.ggrd.comp.d,&tadd,
- FALSE);
- else
- tadd = 0.0;
- }else{ /* cartesian interpolation */
- if(E->X[2][node] <= E->control.slab_theta_bound){
- ggrd_grdtrack_interpolate_xy((double)E->X[1][node],
- (double)E->X[3][node],
- E->control.ggrd.comp.d,&tadd,
- FALSE);
- }else{
- tadd = 0.0;
+ for(hit = slice=0;(!hit) && (slice < E->control.ggrd_slab_slice);slice++){
+ if(E->control.Rsphere) {
+ if(in_slab_slice(E->SX[1][node],slice,E)){
+ /* spherical interpolation */
+ ggrd_grdtrack_interpolate_xy((double)E->SX[2][node] * ONEEIGHTYOVERPI,
+ (double)E->SX[1][node],
+ (E->control.ggrd_ss_grd+slice),&tadd,
+ FALSE);
+ hit = 1;
+ }
+ }else{ /* cartesian interpolation */
+ if(in_slab_slice(E->X[2][node],slice,E)){
+
+ ggrd_grdtrack_interpolate_xy((double)E->X[1][node],
+ (double)E->X[3][node],
+ (E->control.ggrd_ss_grd+slice),&tadd,
+ FALSE);
+ hit = 1 ;
+ }
}
- }
+ } /* end slice loop */
+ if(!hit)
+ tadd = 0;
+ /* end slab slice */
}else{
/* 3-D */
if(E->control.Rsphere) /* spherical interpolation */
@@ -437,7 +462,6 @@
E->control.ggrd.comp.d,&tadd,
FALSE);
}
-
E->C[node] = E->control.ggrd.comp.offset + tadd *
E->control.ggrd.comp.scale;
}
@@ -487,7 +511,56 @@
return;
}
+int in_slab_slice(float coord, int slice, struct All_variables *E)
+{
+ if((slice < 0)||(slice > E->control.ggrd_slab_slice-1))
+ myerror("slab slice out of bounds",E);
+ if(E->control.ggrd_slab_slice < 1)
+ myerror("total slab slice out of bounds",E);
+ if(E->control.ggrd_slab_slice == 1)
+ if(coord <= E->control.ggrd_slab_theta_bound[0])
+ return 1;
+ else
+ return 0;
+ else{
+ if(slice == 0)
+ if(coord <= E->control.ggrd_slab_theta_bound[0])
+ return 1;
+ else
+ return 0;
+ else
+ if((coord <= E->control.ggrd_slab_theta_bound[slice]) && (coord > E->control.ggrd_slab_theta_bound[slice-1]))
+ return 1;
+ else
+ return 0;
+ }
+}
+
+/*
+
+solve a eigenproblem for a symmetric [3][3] matrix (which will not be
+overwritten)
+
+on output, d has the sorted eigenvalues,
+and e the eigenvectors in each column
+
+d[0] > d[1] > d[2]
+
+
+ */
+void ggrd_solve_eigen3x3(double a[3][3],double d[3],
+ double e[3][3],struct All_variables *E)
+{
+ GMT_LONG n=3,nrots;
+ double x[3],b[3],z[3],v[9],a9[9];
+ get_9vec_from_3x3(a9,a);
+ /* the j-th column of V is the eigenvector corresponding to the j-th eigenvalue */
+ if (GMT_jacobi (a9, &n, &n, d, v, b, z, &nrots))
+ myerror("GMT Eigenvalue routine failed to converge in 50 sweeps", E);
+ get_3x3_from_9vec(e,v);
+}
+
#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
/*
@@ -514,12 +587,12 @@
float base[9];
char tfilename[1000];
static ggrd_boolean shift_to_pos_lon = FALSE;
- const int dims=E->mesh.nsd;
+ 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];
-
+ static int init = FALSE;
nox=E->mesh.nox;noy=E->mesh.noy;noz=E->mesh.noz;
elx=E->lmesh.elx;elz=E->lmesh.elz;ely=E->lmesh.ely;
@@ -529,8 +602,12 @@
if(E->viscosity.allow_anisotropic_viscosity == 0)
myerror("ggrd_read_anivisc_from_file: called, but allow_anisotropic_viscosity is FALSE?!",E);
-
- /* isotropic default */
+ if(init)
+ myerror("ggrd_read_anivisc_from_file: called twice",E);
+
+ /*
+ isotropic default
+ */
for(i=E->mesh.levmin;i <= E->mesh.levmax;i++){
nel = E->lmesh.NEL[i];
for(k=1;k <= nel;k++){
@@ -538,6 +615,8 @@
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;
+ E->avmode[i][ind] = (unsigned char)
+ E->viscosity.allow_anisotropic_viscosity;
}
}
}
@@ -563,12 +642,12 @@
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],
+ E->monitor.length_scale*E->viscosity.zbase_layer[E->viscosity.anivisc_layer - 1]/1000,
("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],
+ E->monitor.length_scale/1000*((E->viscosity.anivisc_layer<-1)?(E->viscosity.zbase_layer[-E->viscosity.anivisc_layer - 2]):(0)),
+ E->monitor.length_scale/1000*E->viscosity.zbase_layer[-E->viscosity.anivisc_layer - 1],
("regional"));
/*
@@ -731,8 +810,9 @@
ggrd_grdtrack_free_gstruc(nr_grd);
ggrd_grdtrack_free_gstruc(ntheta_grd);
ggrd_grdtrack_free_gstruc(nphi_grd);
+
+ init = TRUE;
-
}
Modified: mc/3D/CitcomCU/trunk/src/Global_operations.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Global_operations.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Global_operations.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -313,7 +313,29 @@
}
+/* return ||A||^2, where A_i is \int{div(u) d\Omega_i} */
+double global_div_norm2(struct All_variables *E, double *A)
+{
+ int i;
+ double prod, temp;
+ temp = 0.0;
+ prod = 0.0;
+ for (i=1; i<=E->lmesh.npno; i++) {
+ /* L2 norm of div(u) */
+ temp += A[i] * A[i] / E->eco[i].area;
+
+ /* L1 norm */
+ /*temp += fabs(A[m][i]);*/
+ }
+
+ MPI_Allreduce(&temp, &prod, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+
+ return (prod/E->mesh.volume);
+}
+
+
+
double global_vdot(struct All_variables *E, double *A, double *B, int lev)
{
int i, neq;
Modified: mc/3D/CitcomCU/trunk/src/Instructions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Instructions.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Instructions.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -289,6 +289,7 @@
nel = E->lmesh.NEL[i];
nno = E->lmesh.NNO[i];
E->EVI2[i] = (float *) malloc((nel+1)*vpoints[E->mesh.nsd]*sizeof(float));
+ E->avmode[i] = (unsigned char *) malloc((nel+1)*vpoints[E->mesh.nsd]*sizeof(unsigned char));
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));
@@ -306,10 +307,6 @@
}
}
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
@@ -678,7 +675,6 @@
m = E->parallel.me;
- if(E->parallel.me == 0)fprintf(stderr, "ok20\n");
input_string("Problem", E->control.PROBLEM_TYPE, NULL, m);
if(strcmp(E->control.PROBLEM_TYPE, "convection") == 0)
@@ -700,7 +696,6 @@
E->control.CONVECTION = 1;
set_convection_defaults(E);
}
- if(E->parallel.me == 0)fprintf(stderr, "ok21\n");
input_string("Geometry", E->control.GEOMETRY, NULL, m);
if(strcmp(E->control.GEOMETRY, "cart2d") == 0)
@@ -840,8 +835,10 @@
input_string("ggrd_cinit_dfile",E->control.ggrd.comp.dfile,"", m);
input_double("ggrd_cinit_offset",&(E->control.ggrd.comp.offset),"0.0", m);
/* slab slice handling */
- input_boolean("slab_slice",&(E->control.slab_slice),"off", m);
- input_float("slab_theta_bound",&(E->control.slab_theta_bound),"1.0", m);
+ input_int("slab_slice",&(E->control.ggrd_slab_slice),"0", m);
+ if(E->control.ggrd_slab_slice > 3)
+ myerror("too many slab slices",E);
+ input_float_vector("slab_theta_bound",E->control.ggrd_slab_slice,(E->control.ggrd_slab_theta_bound), m);
#endif
@@ -883,7 +880,13 @@
input_float("toptbcval", &(E->control.TBCtopval), "0.0", m);
input_float("bottbcval", &(E->control.TBCbotval), "1.0", m);
+ /* this used to override the VBC settings, I didn't think this
+ should work that way TWB */
input_float("plate_velocity", &(E->control.plate_vel), "0.0", m);
+ if((E->parallel.me == 0) && (fabs(E->control.plate_vel) > 5e-7))
+ fprintf(stderr,"WARNING: plate velocity is overriding VBx \n");
+
+
input_float("plate_age", &(E->control.plate_age), "0.0", m);
input_float("plume_radius", &(E->segment.plume_radius), "0.0", m);
input_float("plume_DT", &(E->segment.plume_DT), "0.0", m);
@@ -1003,11 +1006,8 @@
- 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
===================================================================
--- mc/3D/CitcomCU/trunk/src/Makefile 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Makefile 2010-10-11 04:21:32 UTC (rev 17256)
@@ -164,10 +164,10 @@
-HEADER = element_definitions.h\
- global_defs.h\
- viscosity_descriptions.h\
- advection.h
+HEADER = element_definitions.h \
+ global_defs.h \
+ viscosity_descriptions.h \
+ advection.h prototypes.h
FFILES=#Blas_lapack_interfaces.F
Modified: mc/3D/CitcomCU/trunk/src/Makefile.gzdir
===================================================================
--- mc/3D/CitcomCU/trunk/src/Makefile.gzdir 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Makefile.gzdir 2010-10-11 04:21:32 UTC (rev 17256)
@@ -18,8 +18,8 @@
#COMPRESS=/usr/bin/compress
COMPRESS=/bin/gzip
-GMT_DIR = ../../../GMT4.5.3/
-NETCDF_DIR = ../../../netcdf-3.6.2/
+GMT_DIR = $(GMTHOME)
+NETCDF_DIR = $(NETCDFHOME)
HC_DIR = ../../../hc-svn/
LIB_PATH= -L$(HC_DIR)/objects/x86_64/ -L$(GMT_DIR)/lib/ -L$(NETCDF_DIR)/lib/
#LIB_LIST= -lmpi
@@ -172,10 +172,11 @@
-HEADER = element_definitions.h\
- global_defs.h\
- viscosity_descriptions.h\
- advection.h
+HEADER = element_definitions.h \
+ global_defs.h \
+ viscosity_descriptions.h \
+ advection.h \
+ prototypes.h
FFILES=#Blas_lapack_interfaces.F
Modified: mc/3D/CitcomCU/trunk/src/Makefile.gzdir_ani
===================================================================
--- mc/3D/CitcomCU/trunk/src/Makefile.gzdir_ani 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Makefile.gzdir_ani 2010-10-11 04:21:32 UTC (rev 17256)
@@ -18,16 +18,23 @@
#COMPRESS=/usr/bin/compress
COMPRESS=/bin/gzip
-GMT_DIR = ../../../GMT4.5.3/
-NETCDF_DIR = ../../../netcdf-3.6.2/
+GMT_DIR = $(GMTHOME)
+NETCDF_DIR = $(NETCDFHOME)
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
+
+# use the expokit stuff?
+EXPO_DEFINE = -DCITCOM_USE_EXPOKIT
+EXPO_FILES = expokit/dgpadm.f
+EXPO_LIBS = $(EXPO_OBJS) -llapack -lblas
+
+
+LIB= $(LIB_PATH) $(LIB_LIST) $(EXPO_LIBS) -lm
+
DEFINES = -DUSE_GZDIR -DUSE_GGRD -DCITCOM_ALLOW_ANISOTROPIC_VISC \
- -I$(HC_DIR)/ -Werror \
+ -I$(HC_DIR)/ $(EXPO_DEFINE) \
-I$(GMT_DIR)/include/ -I$(NETCDF_DIR)/include/\
-DPREM_MODEL_FILE=\"$(PWD)/prem/prem.dat\"
@@ -40,9 +47,10 @@
#CC=/usr/lib/cmplrs/cc/gemc_cc
#CC=$(HOME)/usr/local/mpich/bin/mpicc
CC=mpicc
-F77=f77
+F77=mpif77
CPP=
+
CEXT=c
FEXT=F # which implies further action of cpp
OBJEXT=o
@@ -114,7 +122,8 @@
LinuxFLAGS=
LinuxLDFLAGS=
#LinuxOPTIM=-g
-LinuxOPTIM=-O2 $(DEFINES)
+#LinuxOPTIM=-O3 -mtune=core2 -Werror $(DEFINES)
+LinuxOPTIM=-O2 -Werror $(DEFINES)
####################################
#PARAGON
@@ -174,30 +183,32 @@
-HEADER = element_definitions.h\
- global_defs.h\
- viscosity_descriptions.h\
- anisotropic_viscosity.h\
- advection.h
+HEADER = element_definitions.h \
+ global_defs.h \
+ viscosity_descriptions.h \
+ anisotropic_viscosity.h \
+ advection.h \
+ prototypes.h
-FFILES=#Blas_lapack_interfaces.F
-OBJFILES=$(CFILES:.c=.o) $(FFILES:.f=.o)
+FFILES=$(EXPO_FILES) #Blas_lapack_interfaces.F
+OBJFILES=$(CFILES:.c=.o) $(FFILES:.f=.o)
+
####################################################################
# Makefile rules follow
####################################################################
-default: citcom.mpi
+default: citcom.mpi
citcom.mpi: $(OBJFILES) $(HEADER) Makefile
- $(CC) $(OPTIM) $(FLAGS) $(LDFLAGS) -o citcom.mpi \
+ $(F77) $(OPTIM) $(FLAGS) $(LDFLAGS) -o citcom.mpi \
$(OBJFILES) $(FFTLIB) $(LIB)
clean:
- rm -f *.o
+ rm -f $(OBJFILES)
clean-all:
rm -f *.o citcom.mpi
@@ -318,4 +329,5 @@
Ggrd_handling.o: $(HEADER) Ggrd_handling.c
$(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Ggrd_handling.c
-
+expokit/dgpadm.o: expokit/dgpadm.f $(HEADER)
+ $(F77) $(OPTIM) $(FOPTIM) -c expokit/dgpadm.f -o expokit/dgpadm.o
Modified: mc/3D/CitcomCU/trunk/src/Nodal_mesh.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Nodal_mesh.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Nodal_mesh.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -54,7 +54,6 @@
{
//int lev, nodel, i, j, k, ii, d, node;
int lev, nodel, i, j, k, d;
- //float x00, *XX[4], *XG[4], dx[4], dxx[40], dx1, dx2;
float *XX[4], *XG[4], dx[4], dxx[40];
//int n00, nox, noz, noy, fn;
int nox, noz, noy;
Modified: mc/3D/CitcomCU/trunk/src/Output_gzdir.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -76,31 +76,30 @@
int el, els, i, j, k, ii, m, node, fd,snode;
int nox, noz, noy, nfx, nfz, nfy1, nfy2, size1, size2,offset;
static int vtkout = 1;
-
char output_file[1000];
float cvec[3],locx[3];
+ double rtf[4][9];
+ double VV[4][9],lgrad[3][3],isa[3],evel[3];
+ static struct CC Cc;
+ static struct CCX Ccx;
+
static float *SV, *EV;
static int been_here = 0;
float vs;
static int vtk_base_init = 0; /* for spherical */
- static int vtk_pressure_out = 1,
+ const int vtk_pressure_out = 1,
+ vtk_vgm_out = 1,
vtk_viscosity_out = 1;
int vtk_comp_out;
-
-
const int dims = E->mesh.nsd;
const int ends = enodes[dims];
-
const int nno = E->mesh.nno;
+ const int lev = E->mesh.levmax;
+ const int ppts = ppoints[dims];
gzFile gzout;
-
-
-
-
-
if(E->control.composition)
vtk_comp_out = 1;
else
@@ -108,15 +107,12 @@
/* make a directory */
if(E->parallel.me == 0){
- fprintf(stderr,"Output_gzdir: processing output\n");
-
-
-
+ fprintf(stderr,"Output_gzdir: output directory");
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);
+ fprintf(stderr," %s...",output_file);
system(output_file);
- //fprintf(stderr,"making directory: %s\n",output_file);
-
+ fprintf(stderr,"mkdir done\n");
}
/* and wait for the other jobs */
parallel_process_sync();
@@ -144,11 +140,11 @@
gzprintf(gzout, "%6d\n", E->lmesh.nno);
if(!E->control.Rsphere)
for(i = 1; i <= E->lmesh.nno; i++){
- gzprintf(gzout, "%.5e %.5e %.5e\n", E->X[1][i], E->X[2][i], E->X[3][i]);
+ gzprintf(gzout, "%13.6e %13.6e %13.6e\n", E->X[1][i], E->X[2][i], E->X[3][i]);
}
else
for(i = 1; i <= E->lmesh.nno; i++)
- gzprintf(gzout, "%.5e %.5e %.5e\n", E->SX[1][i], E->SX[2][i], E->SX[3][i]);
+ gzprintf(gzout, "%13.6e %13.6e %13.6e\n", E->SX[1][i], E->SX[2][i], E->SX[3][i]);
gzclose(gzout);
/*
surface nodes
@@ -159,9 +155,9 @@
for(i = 1; i <= E->lmesh.nsf; i++){
node = E->surf_node[i];
if(!E->control.Rsphere) /* leave in redundancy for now */
- gzprintf(gzout, "%.5e %.5e %.5e\n", E->X[1][node], E->X[2][node], E->X[3][node]);
+ gzprintf(gzout, "%13.6e %13.6e %13.6e\n", E->X[1][node], E->X[2][node], E->X[3][node]);
else
- gzprintf(gzout, "%.5e %.5e %.5e\n", E->SX[1][node], E->SX[2][node], E->SX[3][node]);
+ gzprintf(gzout, "%13.6e %13.6e %13.6e\n", E->SX[1][node], E->SX[2][node], E->SX[3][node]);
}
gzclose(gzout);
}else{ /* regular I/O */
@@ -175,10 +171,10 @@
fprintf(E->filed[13], "%6d\n", E->lmesh.nno);
if(!E->control.Rsphere)
for(i = 1; i <= E->lmesh.nno; i++)
- fprintf(E->filed[13], "%.5e %.5e %.5e\n", E->X[1][i], E->X[2][i], E->X[3][i]);
+ fprintf(E->filed[13], "%13.6e %13.6e %13.6e\n", E->X[1][i], E->X[2][i], E->X[3][i]);
else
for(i = 1; i <= E->lmesh.nno; i++)
- fprintf(E->filed[13], "%.5e %.5e %.5e\n", E->SX[1][i], E->SX[2][i], E->SX[3][i]);
+ fprintf(E->filed[13], "%13.6e %13.6e %13.6e\n", E->SX[1][i], E->SX[2][i], E->SX[3][i]);
fclose(E->filed[13]);
if((E->parallel.me_loc[3] == E->parallel.nprocz - 1) || /* top */
(E->parallel.me_loc[3] == 0)){ /* bottom */
@@ -191,9 +187,9 @@
for(i = 1; i <= E->lmesh.nsf; i++){
node = E->surf_node[snode];
if(!E->control.Rsphere) /* leave in redundancy for now */
- fprintf(E->filed[13], "%.5e %.5e %.5e\n", E->X[1][node], E->X[2][node], E->X[3][node]);
+ fprintf(E->filed[13], "%13.6e %13.6e %13.6e\n", E->X[1][node], E->X[2][node], E->X[3][node]);
else
- fprintf(E->filed[13], "%.5e %.5e %.5e\n", E->SX[1][node], E->SX[2][node], E->SX[3][node]);
+ fprintf(E->filed[13], "%13.6e %13.6e %13.6e\n", E->SX[1][node], E->SX[2][node], E->SX[3][node]);
}
fclose(E->filed[13]);
}
@@ -257,7 +253,7 @@
E->parallel.me, file_number);
E->filed[10] = safe_fopen(output_file, "w");
/* header line */
- fprintf(E->filed[10], "# %6d %6d %.5e %.5e %.5e %.4e %.4e %.5e %.5e\n",
+ fprintf(E->filed[10], "# %6d %6d %13.6e %13.6e %13.6e %12.4e %12.4e %13.6e %13.6e\n",
E->lmesh.noz, E->advection.timesteps, E->monitor.elapsed_time,
E->slice.Nut, E->slice.Nub, E->data.T_adi0, E->data.T_adi1,
E->monitor.Sigma_interior, E->monitor.Sigma_max);
@@ -266,13 +262,13 @@
*/
if(!E->control.Rsphere){ /* cartesian */
for(i = 1; i <= E->lmesh.noz; i++){
- fprintf(E->filed[10], "%.4e %.4e %.4e %.4e %.4e %.4e %.4e %.4e %.4e\n",
+ fprintf(E->filed[10], "%12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e\n",
E->X[3][i],E->Have.T[i], E->Have.vrms[i], E->Have.Vi[i], E->Have.Rho[i],
E->Have.F[i], E->Have.f[i], E->Have.C[i], E->Have.Tadi[i]);
}
}else{ /* spherical */
for(i = 1; i <= E->lmesh.noz; i++){
- fprintf(E->filed[10], "%.4e %.4e %.4e %.4e %.4e %.4e %.4e %.4e %.4e\n",
+ fprintf(E->filed[10], "%12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e\n",
E->SX[3][i],E->Have.T[i], E->Have.vrms[i], E->Have.Vi[i], E->Have.Rho[i],
E->Have.F[i], E->Have.f[i], E->Have.C[i], E->Have.Tadi[i]);
}
@@ -292,7 +288,7 @@
sprintf(output_file,"%s/%d/t.%d.%d.gz",
E->control.data_file2,file_number,E->parallel.me,file_number);
gzout=safe_gzopen(output_file,"w");
- gzprintf(gzout,"%d %d %.5e\n",
+ gzprintf(gzout,"%d %d %13.6e\n",
file_number,E->lmesh.nno,
E->monitor.elapsed_time);
gzprintf(gzout,"%3d %7d\n",1,E->lmesh.nno); /* for backward compatibility */
@@ -330,8 +326,9 @@
gzprintf(gzout,"%.6e %.6e %.6e\n",cvec[0],cvec[1],cvec[2]);
}
}else{ /* cartesian output in cm/yr */
- for(i=1;i<=E->lmesh.nno;i++)
+ for(i=1;i<=E->lmesh.nno;i++) {
gzprintf(gzout,"%.6e %.6e %.6e\n",E->V[1][i],E->V[2][i],E->V[3][i]);
+ }
}
gzclose(gzout);
if(vtk_pressure_out){
@@ -341,12 +338,99 @@
sprintf(output_file,"%s/%d/p.%d.%d.gz",E->control.data_file2,
file_number, E->parallel.me,file_number);
gzout = safe_gzopen(output_file,"w");
- gzprintf(gzout,"%d %d %.5e\n",file_number,E->lmesh.nno,E->monitor.elapsed_time);
+ gzprintf(gzout,"%d %d %13.6e\n",file_number,E->lmesh.nno,E->monitor.elapsed_time);
gzprintf(gzout,"%3d %7d\n",1,E->lmesh.nno);
for(i=1;i<=E->lmesh.nno;i++)
gzprintf(gzout,"%.6e\n",E->NP[i]);
gzclose(gzout);
}
+ if(vtk_vgm_out){
+ /*
+
+ print velocity gradient matrix and ISA
+
+ */
+ sprintf(output_file,"%s/%d/edot.%d.%d.gz",
+ E->control.data_file2,
+ file_number, E->parallel.me,file_number);
+ gzout = safe_gzopen(output_file,"w");
+ for(i=1;i <= E->lmesh.nel;i++){ /* element loop */
+ if(E->control.Rsphere){ /* need rtf for spherical */
+ get_rtf(E, i, 1, rtf, lev); /* pressure points */
+ if((i-1)%E->lmesh.elz==0)
+ construct_c3x3matrix_el(E,i,&Cc,&Ccx,lev,1);
+ }
+ for(j = 1; j <= ends; j++){ /* velocity at element nodes */
+ node = E->ien[i].node[j];
+ VV[1][j] = E->V[1][node];
+ VV[2][j] = E->V[2][node];
+ VV[3][j] = E->V[3][node];
+ }
+ /* find mean location */
+ locx[0] = locx[1] = locx[2] = 0.0;
+ for(j = 1; j <= ends; j++){
+ node = E->ien[i].node[j];
+ if(E->control.Rsphere){
+ rtp2xyz((float)E->SX[3][node],
+ (float)E->SX[1][node],
+ (float)E->SX[2][node],cvec);
+ locx[0] += cvec[0];
+ locx[1] += cvec[1];
+ locx[2] += cvec[2];
+ }else{
+ locx[0] += E->X[1][node];
+ locx[1] += E->X[2][node];
+ locx[2] += E->X[3][node];
+ }
+ }
+ locx[0] /= ends;locx[1] /= ends;locx[2] /= ends;
+ /* calculate velocity gradient matrix */
+ get_vgm_p(VV,&(E->N),&(E->GNX[lev][i]),&Cc, &Ccx,rtf,
+ E->mesh.nsd,ppts,ends,(E->control.Rsphere),lgrad,
+ evel);
+#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
+ /* calculate the ISA axis or whichever was used to align */
+ switch(E->viscosity.anisotropic_init){
+ case 3:
+ calc_isa_from_vgm(lgrad,evel,i,isa,E,CITCOM_ANIVISC_ALIGN_WITH_VEL);
+ break;
+ case 4:
+ calc_isa_from_vgm(lgrad,evel,i,isa,E,CITCOM_ANIVISC_ALIGN_WITH_ISA);
+ break;
+ case 5:
+ calc_isa_from_vgm(lgrad,evel,i,isa,E,CITCOM_ANIVISC_MIXED_ALIGN);
+ break;
+ default:
+ calc_isa_from_vgm(lgrad,evel,i,isa,E,CITCOM_ANIVISC_ALIGN_WITH_ISA);
+ break;
+ }
+ normalize_vec3d(evel,(evel+1),(evel+2));
+#else
+ isa[0]=isa[1]=isa[2]=evel[0]=evel[1]=evel[2]=0.0;
+#endif
+ /* */
+ if(E->control.Rsphere){
+ /* print velocity gradient matrix, ISA, and normalized velocities */
+ xyz2rtp(locx[0],locx[1],locx[2],cvec);
+ /* r theta phi d[upper right half] (CitcomS spherical system) ISA[3] */
+ gzprintf(gzout,"%11g %11g %11g\t%13.6e %13.6e %13.6e %13.6e %13.6e %13.6e %13.6e %13.6e %13.6e\t%13.6e %13.6e %13.6e\t%13.6e %13.6e %13.6e\n",
+ cvec[0],cvec[1],cvec[2],
+ lgrad[0][0],lgrad[0][1],lgrad[0][2],
+ lgrad[1][0],lgrad[1][1],lgrad[1][2],
+ lgrad[2][0],lgrad[2][1],lgrad[2][2],
+ isa[0],isa[1],isa[2],evel[0],evel[1],evel[2]);
+ }else{
+ gzprintf(gzout,"%11g %11g %11g\t%13.6e %13.6e %13.6e %13.6e %13.6e %13.6e %13.6e %13.6e %13.6e\t%13.6e %13.6e %13.6e\t%13.6e %13.6e %13.6e\n",
+ locx[0],locx[1],locx[2],
+ lgrad[0][0],lgrad[0][1],lgrad[0][2],
+ lgrad[1][0],lgrad[1][1],lgrad[1][2],
+ lgrad[2][0],lgrad[2][1],lgrad[2][2],
+ isa[0],isa[1],isa[2],evel[0],evel[1],evel[2]);
+ }
+ } /* end element loop */
+ gzclose(gzout);
+ }
+
if(vtk_viscosity_out){
/*
@@ -359,7 +443,7 @@
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,"%.4e\n",E->VI[E->mesh.levmax][i]);
+ gzprintf(gzout,"%12.4e\n",E->VI[lev][i]);
}
gzclose(gzout);
#ifdef CITCOM_ALLOW_ANISOTROPIC_VISC
@@ -370,10 +454,10 @@
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]);
+ E->VI2[lev][i],
+ E->VIn1[lev][i],
+ E->VIn2[lev][i],
+ E->VIn3[lev][i]);
}
gzclose(gzout);
@@ -394,7 +478,7 @@
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,"%.4e\n",E->C[i]);
+ gzprintf(gzout,"%12.4e\n",E->C[i]);
gzclose(gzout);
}
/*
@@ -410,60 +494,60 @@
if(E->control.COMPRESS){
sprintf(output_file, "%s/temp.%d.%d.gz", E->control.data_file2, E->parallel.me, file_number);
gzout = safe_gzopen(output_file, "w");
- gzprintf(gzout, "%6d %6d %.5e\n", E->lmesh.nno, E->advection.timesteps,
+ gzprintf(gzout, "%6d %6d %13.6e\n", E->lmesh.nno, E->advection.timesteps,
E->monitor.elapsed_time);
if((file_number % 2* E->control.record_every) == 0)
{
if(E->control.composition == 0)
for(i = 1; i <= E->lmesh.nno; i++)
- // gzprintf(gzout,"%.5e %.4e %.4e %.5e %.5e %.5e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->V[1][i],E->V[2][i],E->V[3][i]);
- // gzprintf(gzout,"%.5e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i]);
- gzprintf(gzout, "%.5e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i]);
+ // gzprintf(gzout,"%13.6e %12.4e %12.4e %13.6e %13.6e %13.6e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->V[1][i],E->V[2][i],E->V[3][i]);
+ // gzprintf(gzout,"%13.6e %12.4e %12.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i]);
+ gzprintf(gzout, "%13.6e %12.4e %12.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i]);
else if(E->control.composition)
for(i = 1; i <= E->lmesh.nno; i++)
- // gzprintf(gzout,"%.5e %.4e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->C[i]);
- gzprintf(gzout, "%.5e %.4e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i], E->C[i]);
+ // gzprintf(gzout,"%13.6e %12.4e %12.4e %12.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->C[i]);
+ gzprintf(gzout, "%13.6e %12.4e %12.4e %12.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i], E->C[i]);
}
else
{
if(E->control.composition == 0)
for(i = 1; i <= E->lmesh.nno; i++)
- // gzprintf(gzout,"%.5e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i]);
- gzprintf(gzout, "%.5e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i]);
+ // gzprintf(gzout,"%13.6e %12.4e %12.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i]);
+ gzprintf(gzout, "%13.6e %12.4e %12.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i]);
else if(E->control.composition)
for(i = 1; i <= E->lmesh.nno; i++)
- // gzprintf(gzout,"%.5e %.4e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->C[i]);
- gzprintf(gzout, "%.5e %.4e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i], E->C[i]);
+ // gzprintf(gzout,"%13.6e %12.4e %12.4e %12.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->C[i]);
+ gzprintf(gzout, "%13.6e %12.4e %12.4e %12.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i], E->C[i]);
}
gzclose(gzout);
}else{
sprintf(output_file, "%s/temp.%d.%d", E->control.data_file2, E->parallel.me, file_number);
E->filed[10] = safe_fopen(output_file, "w");
- fprintf(E->filed[10], "%6d %6d %.5e\n", E->lmesh.nno, E->advection.timesteps,
+ fprintf(E->filed[10], "%6d %6d %13.6e\n", E->lmesh.nno, E->advection.timesteps,
E->monitor.elapsed_time);
if(file_number % (2 * E->control.record_every) == 0)
{
if(E->control.composition == 0)
for(i = 1; i <= E->lmesh.nno; i++)
- // fprintf(E->filed[10],"%.5e %.4e %.4e %.5e %.5e %.5e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->V[1][i],E->V[2][i],E->V[3][i]);
- // fprintf(E->filed[10],"%.5e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i]);
- fprintf(E->filed[10], "%.5e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i]);
+ // fprintf(E->filed[10],"%13.6e %12.4e %12.4e %13.6e %13.6e %13.6e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->V[1][i],E->V[2][i],E->V[3][i]);
+ // fprintf(E->filed[10],"%13.6e %12.4e %12.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i]);
+ fprintf(E->filed[10], "%13.6e %12.4e %12.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i]);
else if(E->control.composition)
for(i = 1; i <= E->lmesh.nno; i++)
- // fprintf(E->filed[10],"%.5e %.4e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->C[i]);
- fprintf(E->filed[10], "%.5e %.4e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i], E->C[i]);
+ // fprintf(E->filed[10],"%13.6e %12.4e %12.4e %12.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->C[i]);
+ fprintf(E->filed[10], "%13.6e %12.4e %12.4e %12.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i], E->C[i]);
}
else
{
if(E->control.composition == 0)
for(i = 1; i <= E->lmesh.nno; i++)
- // fprintf(E->filed[10],"%.5e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i]);
- fprintf(E->filed[10], "%.5e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i]);
+ // fprintf(E->filed[10],"%13.6e %12.4e %12.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i]);
+ fprintf(E->filed[10], "%13.6e %12.4e %12.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i]);
else if(E->control.composition)
for(i = 1; i <= E->lmesh.nno; i++)
- // fprintf(E->filed[10],"%.5e %.4e %.4e %.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->C[i]);
- fprintf(E->filed[10], "%.5e %.4e %.4e %.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i], E->C[i]);
+ // fprintf(E->filed[10],"%13.6e %12.4e %12.4e %12.4e\n",E->T[i],E->heatflux[i],E->heatflux_adv[i],E->C[i]);
+ fprintf(E->filed[10], "%13.6e %12.4e %12.4e %12.4e\n", E->T[i], E->V[3][i], E->heatflux_adv[i], E->C[i]);
}
fclose(E->filed[10]);
}
@@ -485,18 +569,18 @@
sprintf(output_file, "%s/%d/th_t.%d.%d",
E->control.data_file2, file_number, E->parallel.me, file_number);
E->filed[11] = safe_fopen(output_file, "w");
- fprintf(E->filed[11], "%6d %6d %.5e %.5e\n", E->lmesh.nsf, E->advection.timesteps, E->monitor.elapsed_time, E->slice.Nut);
+ fprintf(E->filed[11], "%6d %6d %13.6e %13.6e\n", E->lmesh.nsf, E->advection.timesteps, E->monitor.elapsed_time, E->slice.Nut);
for(i = 1; i <= E->lmesh.nsf; i++)
- fprintf(E->filed[11], "%.5e %.5e %.5e %.5e\n",
+ fprintf(E->filed[11], "%13.6e %13.6e %13.6e %13.6e\n",
E->slice.tpg[i], E->slice.shflux[i], E->Fas410_b[i], E->Fas670_b[i]);
fclose(E->filed[11]);
}else{
sprintf(output_file, "%s/%d/th_t.%d.%d.gz",
E->control.data_file2, file_number, E->parallel.me, file_number);
gzout = safe_gzopen(output_file, "w");
- gzprintf(gzout, "%6d %6d %.5e %.5e\n", E->lmesh.nsf, E->advection.timesteps, E->monitor.elapsed_time, E->slice.Nut);
+ gzprintf(gzout, "%6d %6d %13.6e %13.6e\n", E->lmesh.nsf, E->advection.timesteps, E->monitor.elapsed_time, E->slice.Nut);
for(i = 1; i <= E->lmesh.nsf; i++)
- gzprintf(gzout, "%.5e %.5e %.5e %.5e\n",
+ gzprintf(gzout, "%13.6e %13.6e %13.6e %13.6e\n",
E->slice.tpg[i], E->slice.shflux[i], E->Fas410_b[i], E->Fas670_b[i]);
gzclose(gzout);
}
@@ -507,19 +591,19 @@
sprintf(output_file, "%s/%d/th_b.%d.%d",
E->control.data_file2, file_number, E->parallel.me, file_number);
E->filed[11] = safe_fopen(output_file, "w");
- fprintf(E->filed[11], "%6d %6d %.5e %.5e\n", E->lmesh.nsf, E->advection.timesteps, E->monitor.elapsed_time, E->slice.Nub);
+ fprintf(E->filed[11], "%6d %6d %13.6e %13.6e\n", E->lmesh.nsf, E->advection.timesteps, E->monitor.elapsed_time, E->slice.Nub);
for(i = 1; i <= E->lmesh.nsf; i++)
- fprintf(E->filed[11], "%.5e %.5e %.5e %.5e\n",
+ fprintf(E->filed[11], "%13.6e %13.6e %13.6e %13.6e\n",
E->slice.tpgb[i], E->slice.bhflux[i], E->Fas410_b[i], E->Fas670_b[i]);
fclose(E->filed[11]);
}else{
sprintf(output_file, "%s/%d/th_b.%d.%d.gz",
E->control.data_file2, file_number, E->parallel.me, file_number);
gzout = safe_gzopen(output_file, "w");
- gzprintf(gzout, "%6d %6d %.5e %.5e\n",
+ gzprintf(gzout, "%6d %6d %13.6e %13.6e\n",
E->lmesh.nsf, E->advection.timesteps, E->monitor.elapsed_time, E->slice.Nub);
for(i = 1; i <= E->lmesh.nsf; i++)
- gzprintf(gzout, "%.5e %.5e %.5e %.5e\n", E->slice.tpgb[i], E->slice.bhflux[i], E->Fas410_b[i], E->Fas670_b[i]);
+ gzprintf(gzout, "%13.6e %13.6e %13.6e %13.6e\n", E->slice.tpgb[i], E->slice.bhflux[i], E->Fas410_b[i], E->Fas670_b[i]);
gzclose(gzout);
}
}
@@ -539,7 +623,7 @@
if(E->control.COMPRESS){ /* compressed output */
sprintf(output_file, "%s/%d/traces.%d.gz", E->control.data_file2, file_number,E->parallel.me);
gzout = safe_gzopen(output_file, "w");
- gzprintf(gzout, "%6d %6d %.5e\n", E->advection.markers, E->advection.timesteps, E->monitor.elapsed_time);
+ gzprintf(gzout, "%6d %6d %13.6e\n", E->advection.markers, E->advection.timesteps, E->monitor.elapsed_time);
for(i = 1; i <= E->advection.markers; i++)
gzprintf(gzout, "%g %g %g %d %d\n", E->XMC[1][i], E->XMC[2][i], E->XMC[3][i], E->CElement[i], E->C12[i]);
for(i = 1; i <= E->lmesh.nel; i++)
@@ -548,7 +632,7 @@
}else{
sprintf(output_file, "%s/%d/traces.%d", E->control.data_file2, file_number, E->parallel.me);
E->filed[10] = safe_fopen(output_file, "w");
- fprintf(E->filed[10], "%6d %6d %.5e\n", E->advection.markers, E->advection.timesteps, E->monitor.elapsed_time);
+ fprintf(E->filed[10], "%6d %6d %13.6e\n", E->advection.markers, E->advection.timesteps, E->monitor.elapsed_time);
for(i = 1; i <= E->advection.markers; i++)
fprintf(E->filed[10], "%g %g %g %d %d\n", E->XMC[1][i], E->XMC[2][i], E->XMC[3][i], E->CElement[i], E->C12[i]);
for(i = 1; i <= E->lmesh.nel; i++)
@@ -558,6 +642,8 @@
}
if(E->parallel.me==0)fprintf(stderr,"vel output done\n");
+ parallel_process_sync();
+
return;
}
Modified: mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -85,13 +85,11 @@
//const int i2 = 670;
H = (float *)malloc((E->lmesh.noz + 1) * sizeof(float));
-
for(i = 1; i <= E->lmesh.nno; i++)
{
j = (i - 1) % (E->lmesh.noz) + 1;
E->buoyancy[i] = E->control.Atemp * E->T[i] * E->expansivity[j] - E->control.Acomp * E->C[i];
}
-
if(E->control.Ra_670 != 0.0 || E->control.Ra_410 != 0.0)
{
@@ -103,9 +101,7 @@
}
}
-
remove_horiz_ave(E, E->buoyancy, H, 0);
-
free((void *)H);
return;
}
@@ -929,9 +925,93 @@
void myerror(char *message, struct All_variables *E)
{
E->control.verbose = 1;
+ fprintf(stderr,"node %3i: error: %s\n",E->parallel.me,message);
record(E,message);
- fprintf(stderr,"node %3i: error: %s\n",
- E->parallel.me,message);
parallel_process_termination();
}
+void get_9vec_from_3x3(double *l,double vgm[3][3])
+{
+ l[0] = vgm[0][0];l[1] = vgm[0][1];l[2] = vgm[0][2];
+ l[3] = vgm[1][0];l[4] = vgm[1][1];l[5] = vgm[1][2];
+ l[6] = vgm[2][0];l[7] = vgm[2][1];l[8] = vgm[2][2];
+}
+
+void get_3x3_from_9vec(double l[3][3], double *l9)
+{
+ l[0][0]=l9[0]; l[0][1]=l9[1]; l[0][2]=l9[2];
+ l[1][0]=l9[3]; l[1][1]=l9[4]; l[1][2]=l9[5];
+ l[2][0]=l9[6]; l[2][1]=l9[7]; l[2][2]=l9[8];
+}
+/* symmetric */
+void get_3x3_from_6vec(double l[3][3], double *l6)
+{
+ l[0][0]=l6[0]; l[0][1]=l6[1]; l[0][2]=l6[2];
+ l[1][0]=l6[1]; l[1][1]=l6[3]; l[1][2]=l6[4];
+ l[2][0]=l6[2]; l[2][1]=l6[4]; l[2][2]=l6[5];
+}
+
+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];
+}
+/*
+ C = A * B
+
+ for 3x3 matrix
+
+*/
+void matmul_3x3(double a[3][3],double b[3][3],double c[3][3])
+{
+ int i,j,k;
+ double tmp;
+ for(i=0;i < 3;i++)
+ for(j=0;j < 3;j++){
+ tmp = 0.;
+ for(k=0;k < 3;k++)
+ tmp += a[i][k] * b[k][j];
+ c[i][j] = tmp;
+ }
+}
+
+void assign_to_3x3(double a[3][3],double val)
+{
+ a[0][0]=a[0][1]=a[0][2]=
+ a[1][0]=a[1][1]=a[1][2]=
+ a[2][0]=a[2][1]=a[2][2] = val;
+}
+void remove_trace_3x3(double a[3][3])
+{
+ double trace;
+ trace = (a[0][0]+a[1][1]+a[2][2])/3;
+ a[0][0] -= trace;
+ a[1][1] -= trace;
+ a[2][2] -= trace;
+}
Modified: mc/3D/CitcomCU/trunk/src/Size_does_matter.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Size_does_matter.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Size_does_matter.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -200,7 +200,8 @@
-void get_rtf(struct All_variables *E, int el, int pressure, double rtf[4][9], int lev)
+void get_rtf(struct All_variables *E, int el,
+ int pressure, double rtf[4][9], int lev)
{
//int i, j, k, d, e;
int i, k, d;
@@ -223,9 +224,9 @@
for(i = 1; i <= ends; i++)
x[d] += E->XX[lev][d][E->IEN[lev][el].node[i]] * E->N.vpt[GNVINDEX(i, k)];
- rtf[3][k] = 1.0 / sqrt(x[1] * x[1] + x[2] * x[2] + x[3] * x[3]);
- rtf[1][k] = acos(x[3] * rtf[3][k]);
- rtf[2][k] = myatan(x[2], x[1]);
+ rtf[3][k] = 1.0 / sqrt(x[1] * x[1] + x[2] * x[2] + x[3] * x[3]); /* 1/r */
+ rtf[1][k] = acos(x[3] * rtf[3][k]); /* theta */
+ rtf[2][k] = myatan(x[2], x[1]); /* phi */
}
}
else
Modified: mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Stokes_flow_Incomp.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -81,17 +81,220 @@
been_here = 1;
v_from_vector(E, E->V, E->U);
- /* p_to_nodes(E,E->P,E->NP,E->mesh.levmax); */
+#ifdef USE_GZDIR
+ if(E->control.gzdir)
+ p_to_nodes(E,E->P,E->NP,E->mesh.levmax);
+#endif
return;
}
-/* ========================================================================== */
+/* ==========================================================================
-float solve_Ahat_p_fhat(struct All_variables *E, double *V, double *P, double *F, double imp, int *steps_max)
+adjusted CitcomS style (doesn't quite work yet)
+
+
+*/
+
+float solve_Ahat_p_fhat_new(struct All_variables *E, double *V, double *P,
+ double *FF, double imp, int *steps_max)
{
+ int i, j, count, convergent, valid, problems, lev, npno, neq;
+ int gnpno, gneq;
+
+ static int been_here = 0;
+ static double *r1, *r2, *z1, *s1, *s2, *u1;
+ double *F;
+ double *shuffle;
+ double alpha, delta, r0dotr0, r1dotz1, r0dotz0;
+ double residual, initial_residual, res_magnitude, v_res;
+ char message[500];
+ double time0, time;
+ static double timea;
+ float dpressure, dvelocity;
+
+ npno = E->lmesh.npno;
+ neq = E->lmesh.neq;
+ lev = E->mesh.levmax;
+ gnpno = E->mesh.npno;
+ gneq = E->mesh.neq;
+
+ if(been_here == 0) {
+ r1 = (double *)malloc((npno + 1) * sizeof(double));
+ r2 = (double *)malloc((npno + 1) * sizeof(double));
+ z1 = (double *)malloc((npno + 1) * sizeof(double));
+ s1 = (double *)malloc((npno + 1) * sizeof(double));
+ s2 = (double *)malloc((npno + 1) * sizeof(double));
+ u1 = (double *)malloc((neq + 2) * sizeof(double));
+ F = (double *)malloc((neq+2)*sizeof(double));
+
+ timea = CPU_time0();
+ }
+ been_here++;
+
+ problems = 0;
+ time0 = time = CPU_time0();
+
+ /* copy the F vector */
+ for(j=0;j<neq;j++)
+ F[j] = FF[j];
+
+ v_res = sqrt(global_vdot(E, E->F, E->F, lev));
+
+
+ initial_vel_residual(E, V, P, F, u1,v_res);
+
+ assemble_div_u(E, V, r1, lev);
+
+ E->monitor.incompressibility = sqrt(global_div_norm2(E, r1)
+ / (1e-32 + E->monitor.vdotv));
+
+ dvelocity = 1.0;
+ dpressure = 1.0;
+ convergent = 0;
+ r0dotz0 = 0;
+
+ if(E->parallel.me == 0)
+ fprintf(stderr, "initial residue of momentum equation %g %d inc %g\n",
+ v_res, gneq,E->monitor.incompressibility );
+
+ count = 0;
+ E->monitor.vdotv = global_vdot(E, V, V, lev);
+ E->monitor.pdotp = global_pdot(E, P, P, lev);
+ E->monitor.vdotv = sqrt(E->monitor.vdotv/gneq);
+ E->monitor.pdotp = sqrt(E->monitor.pdotp/gnpno);
+
+ generate_log_message(count,time0,timea,dvelocity, dpressure,E);
+ residual = 0;
+ while((count < *steps_max) && (dpressure >= imp || dvelocity >= imp)){
+
+ for(j = 1; j <= npno; j++)
+ z1[j] = E->BPI[lev][j] * r1[j];
+
+ r1dotz1 = global_pdot(E, r1, z1, lev);
+ assert(r1dotz1 != 0.0 /* Division by zero in head of incompressibility iteration */);
+
+ if((count == 0))
+ for(j = 1; j <= npno; j++)
+ s2[j] = z1[j];
+ else{
+ delta = r1dotz1 / r0dotr0;
+ for(j = 1; j <= npno; j++)
+ s2[j] = z1[j] + delta * s1[j];
+ }
+
+ assemble_grad_p(E, s2, F, lev);
+
+ valid = solve_del2_u(E, u1, F, imp * v_res, lev);
+ strip_bcs_from_residual(E, u1, lev);
+
+ assemble_div_u(E, u1, F, lev);
+
+ alpha = r1dotz1 / global_pdot(E, s2, F, lev);
+
+ /* r2 = r1 - alpha * div(u1) */
+ for(j=1; j<=npno; j++)
+ r2[j] = r1[j] - alpha * F[j];
+
+ /* P = P + alpha * s2 */
+ for(j=1; j<=npno; j++)
+ P[j] += alpha * s2[j];
+
+ /* V = V - alpha * u1 */
+ for(j=0; j<neq; j++)
+ V[j] -= alpha * u1[j];
+
+ E->monitor.vdotv = global_vdot(E, V, V, lev);
+ E->monitor.pdotp = global_pdot(E, P, P, lev);
+
+ assemble_div_u(E, V, z1, lev);
+ E->monitor.incompressibility = sqrt(global_div_norm2(E, z1)
+ / (1e-32 + E->monitor.vdotv));
+
+ 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++;
+
+ generate_log_message(count,time0,timea,dvelocity, dpressure,E);
+
+ shuffle = s1;
+ s1 = s2;
+ s2 = shuffle;
+
+ shuffle = r1;
+ r1 = r2;
+ r2 = shuffle;
+
+ /* shift <r0, z0> = <r1, z1> */
+ r0dotz0 = r1dotz1;
+
+
+ } /* end loop for conjugate gradient */
+ assemble_div_u(E, V, z1, lev);
+ if(problems){
+ fprintf(E->fp, "Convergence of velocity solver may affect continuity\n");
+ fprintf(E->fp, "Consider running with the `see_convergence=on' option\n");
+ fprintf(E->fp, "To evaluate the performance of the current relaxation parameters\n");
+ fflush(E->fp);
+ }
+
+ if(E->control.print_convergence && E->parallel.me == 0){
+ fprintf(E->fp, "after (%03d) pressure loops and %g sec for step %d\n", count, CPU_time0() - timea, E->monitor.solution_cycles);
+ fprintf(stderr, "after (%03d) pressure loops and %g sec for step %d\n", count, CPU_time0() - timea, E->monitor.solution_cycles);
+ fflush(E->fp);
+ }
+
+
+ *steps_max = count;
+
+ return (residual);
+}
+
+void initial_vel_residual(struct All_variables *E,
+ double *V, double *P, double *F,double *u1,
+ double acc)
+{
+ int neq = E->lmesh.neq;
+ int lev = E->mesh.levmax;
+ int i, valid;
+
+ /* F = F - grad(P) - K*V */
+ assemble_grad_p(E, P, u1, lev);
+ for(i=0; i<neq; i++)
+ F[i] -= u1[i];
+
+ assemble_del2_u(E, V, u1, lev, 1);
+ for(i=0; i<neq; i++)
+ F[i] -= u1[i];
+
+ strip_bcs_from_residual(E, F, lev);
+
+ /* solve K*u1 = F for u1 */
+ valid = solve_del2_u(E, u1, F, acc, lev);
+
+ if(!valid && (E->parallel.me==0)) {
+ fputs("Warning: solver not converging! 0\n", stderr);
+ fputs("Warning: solver not converging! 0\n", E->fp);
+ }
+ strip_bcs_from_residual(E, u1, lev);
+
+ /* V = V + u1 */
+ for(i=0; i < neq; i++)
+ V[i] += u1[i];
+
+ return;
+}
+
+float solve_Ahat_p_fhat(struct All_variables *E, double *V, double *P,
+ double *F, double imp, int *steps_max)
+{
//int i, j, k, ii, count, convergent, valid, problems, lev, lev_low, npno, neq, steps;
int i, j, count, convergent, valid, problems, lev, npno, neq;
int gnpno, gneq;
@@ -119,7 +322,7 @@
gnpno = E->mesh.npno;
gneq = E->mesh.neq;
- if(been_here == 0)
+ if(been_here == 0) /* been here gets incremented further down */
{
r0 = (double *)malloc((npno + 1) * sizeof(double));
r1 = (double *)malloc((npno + 1) * sizeof(double));
@@ -130,7 +333,7 @@
s2 = (double *)malloc((npno + 1) * sizeof(double));
Ah = (double *)malloc((neq + 1) * sizeof(double));
u1 = (double *)malloc((neq + 2) * sizeof(double));
- been_here = 1;
+
timea = CPU_time0();
}
@@ -156,7 +359,7 @@
strip_bcs_from_residual(E, Ah, lev);
- valid = solve_del2_u(E, u1, Ah, imp * v_res, E->mesh.levmax, 0);
+ valid = solve_del2_u(E, u1, Ah, imp * v_res, E->mesh.levmax);
strip_bcs_from_residual(E, u1, lev);
for(i = 0; i < neq; i++)
@@ -209,7 +412,7 @@
assemble_grad_p(E, s2, Ah, lev);
- valid = solve_del2_u(E, u1, Ah, imp * v_res, lev, 1);
+ valid = solve_del2_u(E, u1, Ah, imp * v_res, lev);
strip_bcs_from_residual(E, u1, lev);
assemble_div_u(E, u1, Ah, lev);
@@ -293,6 +496,8 @@
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;
@@ -333,15 +538,15 @@
for(node = 1; node <= nno; node++)
{
- V[1][node] = F[E->id[node].doff[1]];
- V[2][node] = F[E->id[node].doff[2]];
- V[3][node] = F[E->id[node].doff[3]];
- if(E->node[node] & VBX)
- V[1][node] = E->VB[1][node];
- if(E->node[node] & VBZ)
- V[3][node] = E->VB[3][node];
- if(E->node[node] & VBY)
- V[2][node] = E->VB[2][node];
+ V[1][node] = F[E->id[node].doff[1]];
+ V[2][node] = F[E->id[node].doff[2]];
+ V[3][node] = F[E->id[node].doff[3]];
+ if(E->node[node] & VBX)
+ V[1][node] = E->VB[1][node];
+ if(E->node[node] & VBY)
+ V[2][node] = E->VB[2][node];
+ if(E->node[node] & VBZ)
+ V[3][node] = E->VB[3][node];
}
return;
}
Modified: mc/3D/CitcomCU/trunk/src/Topo_gravity.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Topo_gravity.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Topo_gravity.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -354,7 +354,7 @@
Sxz = 0.0;
Szy = 0.0;
if(E->control.Rsphere)
- get_rtf(E, e, 0, rtf, lev);
+ get_rtf(E, e, 0, rtf, lev); /* vpts */
for(j = 1; j <= vpts; j++)
{
Modified: mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2010-10-11 04:21:32 UTC (rev 17256)
@@ -57,7 +57,7 @@
{
int i, l;
float temp;
- int m = E->parallel.me ;
+ int m = E->parallel.me ;
/* default values .... */
E->viscosity.update_allowed = 0;
E->viscosity.SDEPV = E->viscosity.TDEPV = E->viscosity.BDEPV =
@@ -67,7 +67,7 @@
input_boolean("plasticity_dimensional",&(E->viscosity.plasticity_dimensional),"on",m);
- for(i = 0; i < 40; i++)
+ for(i = 0; i < CITCOM_CU_VISC_MAXLAYER; i++)
{
E->viscosity.N0[i] = 1.0;
E->viscosity.T[i] = 0.0;
@@ -98,7 +98,9 @@
/* comp dep visc */
E->viscosity.pre_comp[i] = 1.0;
-
+
+
+ E->viscosity.zbase_layer[i] = -999; /* not assigned by default */
} /* end layer loop */
/* read in information */
@@ -110,8 +112,6 @@
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);
@@ -127,28 +127,20 @@
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)) {
-
+ /* no z_layer input found */
if(E->viscosity.num_mat != 4)
- myerror("error: either use z_layer for non dim layer depths, or set num_mat to four",E);
-
+ myerror("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;
+ if(E->control.Rsphere)
+ E->viscosity.zbase_layer[3] = 0.55;
+ else
+ E->viscosity.zbase_layer[3] = 0.;
}
-
-
-
-
-
-
-
-
-
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);
@@ -194,7 +186,15 @@
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_int("anisotropic_init",&(E->viscosity.anisotropic_init),"0",m); /* 0: isotropic
+ 1: random
+ 2: read in director orientation
+ and log10(eta_s/eta)
+ 3: align with velocity
+ 4: align with ISA
+ 5: align mixed depending on deformation state
+
+ */
input_string("anisotropic_init_dir",(E->viscosity.anisotropic_init_dir),"",m); /* directory
for
ggrd
@@ -209,12 +209,20 @@
from
isotropic
solution? */
+ if(!E->viscosity.anivisc_start_from_iso)
+ if(E->viscosity.anisotropic_init == 3){
+ if(E->parallel.me == 0)fprintf(stderr,"WARNING: overriding anisotropic first step for ani init mode\n");
+ E->viscosity.anivisc_start_from_iso = TRUE;
+ }
+ /* ratio between weak and strong direction */
+ input_double("ani_vis2_factor",&(E->viscosity.ani_vis2_factor),"0.1",m);
+
+
}
#endif
-
/* composition factors */
input_float_vector("pre_comp",2,(E->viscosity.pre_comp),m);
@@ -689,127 +697,279 @@
}
-void strain_rate_2_inv(struct All_variables *E, float *EEDOT, int SQRT)
+void strain_rate_2_inv(struct All_variables *E, float *EEDOT,
+ int SQRT)
{
- double edot[4][4], dudx[4][4], rtf[4][9];
- float VV[4][9], Vxyz[9][9];
+ double edot[4][4], dudx[4][4], rtf[4][9];
+ float VV[4][9], Vxyz[9][9];
+
+ //int e, i, j, p, q, n, nel, k;
+ int e, i, j, p, q, n, nel;
+
+ const int dims = E->mesh.nsd;
+ const int ends = enodes[dims];
+ const int lev = E->mesh.levmax;
+ const int ppts = ppoints[dims];
+
+ nel = E->lmesh.nel;
+
+
+ if(E->control.Rsphere){
+
+ for(e = 1; e <= nel; e++){
- //int e, i, j, p, q, n, nel, k;
- int e, i, j, p, q, n, nel;
+ get_rtf(E, e, 1, rtf, lev); /* pressure */
+
+ for(i = 1; i <= ends; i++){
+ n = E->ien[e].node[i];
+ VV[1][i] = E->V[1][n];
+ VV[2][i] = E->V[2][n];
+ VV[3][i] = E->V[3][n];
+ }
+
+ for(j = 1; j <= ppts; j++){
+ Vxyz[1][j] = 0.0;
+ Vxyz[2][j] = 0.0;
+ Vxyz[3][j] = 0.0;
+ Vxyz[4][j] = 0.0;
+ Vxyz[5][j] = 0.0;
+ Vxyz[6][j] = 0.0;
+ }
- const int dims = E->mesh.nsd;
- const int ends = enodes[dims];
- const int lev = E->mesh.levmax;
- //const int nno = E->lmesh.nno;
- //const int vpts = vpoints[dims];
- const int ppts = ppoints[dims];
- //const int sphere_key = 1;
+ for(j = 1; j <= ppts; j++){ /* only makes "sense" for ppts = 1 */
+ for(i = 1; i <= ends; i++){
+ Vxyz[1][j] += (VV[1][i] * E->gNX[e].ppt[GNPXINDEX(0, i, j)] + VV[3][i] * E->N.ppt[GNPINDEX(i, j)]) * rtf[3][j]; /* tt */
+ Vxyz[2][j] += ((VV[2][i] * E->gNX[e].ppt[GNPXINDEX(1, i, j)] + VV[1][i] * E->N.ppt[GNPINDEX(i, j)] * cos(rtf[1][j])) / sin(rtf[1][j]) + VV[3][i] * E->N.ppt[GNPINDEX(i, j)]) * rtf[3][j]; /* pp */
+ Vxyz[3][j] += VV[3][i] * E->gNX[e].ppt[GNPXINDEX(2, i, j)]; /* rr */
+
+ Vxyz[4][j] += ((VV[1][i] * E->gNX[e].ppt[GNPXINDEX(1, i, j)] - VV[2][i] * E->N.ppt[GNPINDEX(i, j)] * cos(rtf[1][j])) / sin(rtf[1][j]) + VV[2][i] * E->gNX[e].ppt[GNPXINDEX(0, i, j)]) * rtf[3][j]; /* tp */
+ Vxyz[5][j] += VV[1][i] * E->gNX[e].ppt[GNPXINDEX(2, i, j)] + rtf[3][j] * (VV[3][i] * E->gNX[e].ppt[GNPXINDEX(0, i, j)] - VV[1][i] * E->N.ppt[GNPINDEX(i, j)]); /* tr */
+ Vxyz[6][j] += VV[2][i] * E->gNX[e].ppt[GNPXINDEX(2, i, j)] + rtf[3][j] * (VV[3][i] * E->gNX[e].ppt[GNPXINDEX(1, i, j)] / sin(rtf[1][j]) - VV[2][i] * E->N.ppt[GNPINDEX(i, j)]); /* pr */
+ }
+ edot[1][1] = 2.0 * Vxyz[1][j]; /* this should be a summation */
+ edot[2][2] = 2.0 * Vxyz[2][j];
+ edot[3][3] = 2.0 * Vxyz[3][j];
+ edot[1][2] = Vxyz[4][j];
+ edot[1][3] = Vxyz[5][j];
+ edot[2][3] = Vxyz[6][j];
+ }
+ /* /\* eps: rr, rt, rp, tt, tp, pp *\/ */
+ /* if(e < 5)fprintf(stderr,"1: %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e\n",edot[3][3]/2,edot[1][3]/2,edot[2][3]/2,edot[1][1]/2,edot[1][2]/2,edot[2][2]/2); */
+ EEDOT[e] = edot[1][1] * edot[1][1] + edot[1][2] * edot[1][2] * 2.0 + edot[2][2] * edot[2][2] + edot[2][3] * edot[2][3] * 2.0 + edot[3][3] * edot[3][3] + edot[1][3] * edot[1][3] * 2.0;
+
+ }
- nel = E->lmesh.nel;
+ }else if(E->control.CART3D){
+ for(e = 1; e <= nel; e++){
+ for(i = 1; i <= ends; i++){
+ n = E->ien[e].node[i];
+ VV[1][i] = E->V[1][n];
+ VV[2][i] = E->V[2][n];
+ if(dims == 3)
+ VV[3][i] = E->V[3][n];
+ }
+ for(p = 1; p <= dims; p++)
+ for(q = 1; q <= dims; q++)
+ dudx[p][q] = 0.0;
+
+ for(i = 1; i <= ends; i++)
+ for(p = 1; p <= dims; p++)
+ for(q = 1; q <= dims; q++)
+ dudx[p][q] += VV[p][i] * E->gNX[e].ppt[GNPXINDEX(q - 1, i, 1)];
+
+ for(p = 1; p <= dims; p++)
+ for(q = 1; q <= dims; q++)
+ edot[p][q] = dudx[p][q] + dudx[q][p];
- if(E->control.Rsphere)
- {
+ /* if(e < 5)fprintf(stderr,"1: %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e\n",edot[1][1]/2,edot[1][2]/2,edot[1][3]/2,edot[2][2]/2,edot[2][3]/2,edot[3][3]/2); */
- for(e = 1; e <= nel; e++)
- {
+
+ if(dims == 2)
+ EEDOT[e] = edot[1][1] * edot[1][1] + edot[2][2] * edot[2][2] + edot[1][2] * edot[1][2] * 2.0;
+
+ else if(dims == 3)
+ EEDOT[e] = edot[1][1] * edot[1][1] + edot[1][2] * edot[1][2] * 2.0 + edot[2][2] * edot[2][2] + edot[2][3] * edot[2][3] * 2.0 + edot[3][3] * edot[3][3] + edot[1][3] * edot[1][3] * 2.0;
+
+ }
+ }
+
- get_rtf(E, e, 2, rtf, lev);
+
+ if(SQRT)
+ for(e = 1; e <= nel; e++)
+ EEDOT[e] = sqrt(0.5 * EEDOT[e]);
+ else
+ for(e = 1; e <= nel; e++)
+ EEDOT[e] *= 0.5;
+
+ return;
+}
- for(i = 1; i <= ends; i++)
- {
- n = E->ien[e].node[i];
- VV[1][i] = E->V[1][n];
- VV[2][i] = E->V[2][n];
- VV[3][i] = E->V[3][n];
- }
+/* compute second invariant from a strain-rate tensor in 0,...2 format
- for(j = 1; j <= ppts; j++)
- {
- Vxyz[1][j] = 0.0;
- Vxyz[2][j] = 0.0;
- Vxyz[3][j] = 0.0;
- Vxyz[4][j] = 0.0;
- Vxyz[5][j] = 0.0;
- Vxyz[6][j] = 0.0;
- }
+ */
+double second_invariant_from_3x3(double e[3][3])
+{
+ return(sqrt(0.5*
+ (e[0][0] * e[0][0] +
+ e[0][1] * e[0][1] * 2.0 +
+ e[1][1] * e[1][1] +
+ e[1][2] * e[1][2] * 2.0 +
+ e[2][2] * e[2][2] +
+ e[0][2] * e[0][2] * 2.0)));
+}
- for(j = 1; j <= ppts; j++)
- {
- for(i = 1; i <= ends; i++)
- {
- Vxyz[1][j] += (VV[1][i] * E->gNX[e].ppt[GNPXINDEX(0, i, j)] + VV[3][i] * E->N.ppt[GNPINDEX(i, j)]) * rtf[3][j];
- Vxyz[2][j] += ((VV[2][i] * E->gNX[e].ppt[GNPXINDEX(1, i, j)] + VV[1][i] * E->N.ppt[GNPINDEX(i, j)] * cos(rtf[1][j])) / sin(rtf[1][j]) + VV[3][i] * E->N.ppt[GNPINDEX(i, j)]) * rtf[3][j];
- Vxyz[3][j] += VV[3][i] * E->gNX[e].ppt[GNPXINDEX(2, i, j)];
+/*
- Vxyz[4][j] += ((VV[1][i] * E->gNX[e].ppt[GNPXINDEX(1, i, j)] - VV[2][i] * E->N.ppt[GNPINDEX(i, j)] * cos(rtf[1][j])) / sin(rtf[1][j]) + VV[2][i] * E->gNX[e].ppt[GNPXINDEX(0, i, j)]) * rtf[3][j];
- Vxyz[5][j] += VV[1][i] * E->gNX[e].ppt[GNPXINDEX(2, i, j)] + rtf[3][j] * (VV[3][i] * E->gNX[e].ppt[GNPXINDEX(0, i, j)] - VV[1][i] * E->N.ppt[GNPINDEX(i, j)]);
- Vxyz[6][j] += VV[2][i] * E->gNX[e].ppt[GNPXINDEX(2, i, j)] + rtf[3][j] * (VV[3][i] * E->gNX[e].ppt[GNPXINDEX(1, i, j)] / sin(rtf[1][j]) - VV[2][i] * E->N.ppt[GNPINDEX(i, j)]);
- }
- edot[1][1] = 2.0 * Vxyz[1][j];
- edot[2][2] = 2.0 * Vxyz[2][j];
- edot[3][3] = 2.0 * Vxyz[3][j];
- edot[1][2] = Vxyz[4][j];
- edot[1][3] = Vxyz[5][j];
- edot[2][3] = Vxyz[6][j];
- }
+ calculate velocity gradient matrix for all elements
- EEDOT[e] = edot[1][1] * edot[1][1] + edot[1][2] * edot[1][2] * 2.0 + edot[2][2] * edot[2][2] + edot[2][3] * edot[2][3] * 2.0 + edot[3][3] * edot[3][3] + edot[1][3] * edot[1][3] * 2.0;
+ l is stored as (e-1)*9
+
+*/
+void calc_vgm_matrix(struct All_variables *E, double *l,double *evel)
+{
+ double rtf[4][9];
+ double VV[4][9],vgm[3][3];
+ int e,j,i,n,p,q,eoff;
+ double d[6];
+ static struct CC Cc;
+ static struct CCX Ccx;
+ const int dims = E->mesh.nsd;
+ const int ppts = ppoints[dims];
+ const int ends = enodes[dims];
+ const int lev = E->mesh.levmax;
+ const int nel = E->lmesh.nel;
- }
+ for(eoff=0,e=1; e <= nel; e++, eoff += 9) {
+ if(E->control.Rsphere){ /* need rtf for spherical */
+ get_rtf(E, e, 1, rtf, lev); /* pressure points */
+ if((e-1)%E->lmesh.elz==0)
+ construct_c3x3matrix_el(E,e,&Cc,&Ccx,lev,1);
+ }
+ for(i = 1; i <= ends; i++){ /* velocity at element nodes */
+ n = E->ien[e].node[i];
+ VV[1][i] = E->V[1][n];
+ VV[2][i] = E->V[2][n];
+ VV[3][i] = E->V[3][n];
+ }
+ get_vgm_p(VV,&(E->N),&(E->GNX[lev][e]),&Cc, &Ccx,rtf,
+ E->mesh.nsd,ppts,ends,(E->control.Rsphere),vgm,evel);
+ get_9vec_from_3x3((l+eoff),vgm);
+ }
+}
+/*
- }
+ given a 3x3 velocity gradient matrix l, compute a d[3][3]
+ strain-rate matrix
- else if(E->control.CART3D)
- {
+*/
- for(e = 1; e <= nel; e++)
- {
+void calc_strain_from_vgm(double l[3][3], double d[3][3])
+{
+ int i,j;
+ for(i=0;i < 3;i++)
+ for(j=0;j < 3;j++)
+ d[i][j] = 0.5 * (l[i][j] + l[j][i]);
+}
+void calc_strain_from_vgm9(double *l9, double d[3][3])
+{
+ double l[3][3];
+ get_3x3_from_9vec(l, l9);
+ calc_strain_from_vgm(l, d);
+}
+/*
- for(i = 1; i <= ends; i++)
- {
- n = E->ien[e].node[i];
- VV[1][i] = E->V[1][n];
- VV[2][i] = E->V[2][n];
- if(dims == 3)
- VV[3][i] = E->V[3][n];
- }
+ given a 3x3 velocity gradient matrix l, compute a rotation matrix
- for(p = 1; p <= dims; p++)
- for(q = 1; q <= dims; q++)
- dudx[p][q] = 0.0;
+*/
- for(i = 1; i <= ends; i++)
- for(p = 1; p <= dims; p++)
- for(q = 1; q <= dims; q++)
- dudx[p][q] += VV[p][i] * E->gNX[e].ppt[GNPXINDEX(q - 1, i, 1)];
+void calc_rot_from_vgm(double l[3][3], double r[3][3])
+{
+ int i,j;
+ for(i=0;i < 3;i++)
+ for(j=0;j < 3;j++)
+ r[i][j] = 0.5 * (l[i][j] - l[j][i]);
+}
- for(p = 1; p <= dims; p++)
- for(q = 1; q <= dims; q++)
- edot[p][q] = dudx[p][q] + dudx[q][p];
+/*
- if(dims == 2)
- EEDOT[e] = edot[1][1] * edot[1][1] + edot[2][2] * edot[2][2] + edot[1][2] * edot[1][2] * 2.0;
+ different version of strain-rate computation, obtain strain-rate
+ matrix at all elements, storing it in edot as (e-1)*6 upper
+ triangle format
- else if(dims == 3)
- EEDOT[e] = edot[1][1] * edot[1][1] + edot[1][2] * edot[1][2] * 2.0 + edot[2][2] * edot[2][2] + edot[2][3] * edot[2][3] * 2.0 + edot[3][3] * edot[3][3] + edot[1][3] * edot[1][3] * 2.0;
- }
- }
+*/
+void calc_strain_rate_matrix(struct All_variables *E,
+ double *edot)
+{
+ double rtf[4][9];
+ double VV[4][9], Vxyz[7][9];
+ int e,j,i,n,p,q,eoff;
+ static struct CC Cc;
+ static struct CCX Ccx;
+ double ba[9][4][9][7];
+ const int dims = E->mesh.nsd;
+ const int ppts = ppoints[dims];
+ const int ends = enodes[dims];
+ const int lev = E->mesh.levmax;
+ const int nel = E->lmesh.nel;
+ //fprintf(stderr,"\n");
+ for(eoff=0,e=1; e <= nel; e++,eoff+=6) {
+ if(E->control.Rsphere){ /* need rtf for spherical */
+ get_rtf(E, e, 1, rtf, lev); /* pressure */
+ if((e-1)%E->lmesh.elz==0)
+ construct_c3x3matrix_el(E,e,&Cc,&Ccx,lev,1);
+ }
+ for(i = 1; i <= ends; i++){ /* velocity at element nodes */
+ n = E->ien[e].node[i];
+ VV[1][i] = E->V[1][n];
+ VV[2][i] = E->V[2][n];
+ VV[3][i] = E->V[3][n];
+ }
+ for(j = 1; j <= ppts; j++){
+ Vxyz[1][j] = 0.0;
+ Vxyz[2][j] = 0.0;
+ Vxyz[3][j] = 0.0;
+ Vxyz[4][j] = 0.0;
+ Vxyz[5][j] = 0.0;
+ Vxyz[6][j] = 0.0;
+ }
+ get_ba_p(&(E->N),&(E->GNX[lev][e]),&Cc, &Ccx,rtf,
+ E->mesh.nsd,ppts,ends,(E->control.Rsphere),ba);
+ for(j=1;j <= ppts;j++)
+ for(p=1;p <= 6;p++)
+ for(i=1;i <= ends;i++)
+ for(q=1;q <= dims;q++) {
+ Vxyz[p][j] += ba[i][q][j][p] * VV[q][i];
+ }
+ edot[eoff+0] = edot[eoff+1] = edot[eoff+2] =
+ edot[eoff+3] = edot[eoff+4] = edot[eoff+5] = 0.0;
+ for(j=1; j <= ppts; j++) {
+ edot[eoff+0] += Vxyz[1][j]; /* e_xx = e_tt */
+ edot[eoff+1] += Vxyz[4][j]/2; /* e_xy = e_tp */
+ edot[eoff+2] += Vxyz[5][j]/2; /* e_xz = e_tr */
+ edot[eoff+3] += Vxyz[2][j]; /* e_yy = e_pp */
+ edot[eoff+4] += Vxyz[6][j]/2; /* e_yz = e_pr */
+ edot[eoff+5] += Vxyz[3][j]; /* e_zz = e_rr */
+ }
+ if(ppts != 1){
+ for(i=0;i<6;i++)
+ edot[eoff+i] /= (float)ppts;
+ }
- if(SQRT)
- for(e = 1; e <= nel; e++)
- EEDOT[e] = sqrt(0.5 * EEDOT[e]);
- else
- for(e = 1; e <= nel; e++)
- EEDOT[e] *= 0.5;
-
- return;
+ /* if(e < 5){ */
+ /* if(E->control.Rsphere){ */
+ /* /\* rr rt rp tt tp pp *\/ */
+ /* fprintf(stderr,"2: %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e \n",edot[eoff+5],edot[eoff+2],edot[eoff+4],edot[eoff+0],edot[eoff+1],edot[eoff+3]); */
+ /* }else{ */
+ /* fprintf(stderr,"2: %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e \n",edot[eoff+0],edot[eoff+1],edot[eoff+2],edot[eoff+3],edot[eoff+4],edot[eoff+5]); */
+ /* } */
+ /* } */
+ }
}
-
-
int layers(struct All_variables *E, float x3)
{
int llayers = 0;
@@ -979,7 +1139,7 @@
eta_old2 = eta_old * eta_old;
eta_new = (tau2 * eta_old)/(tau2 + 2.0 * eta_old2 * eedot[i]);
EEta[ (i-1)*vpts + jj ] = ettnew;
- //if(E->parallel.me==0)fprintf(stderr,"tau: %11g eII: %11g eta_old: %11g eta_new: %11g\n",tau, eedot[i],eta_old,eta_new);
+ //if(E->parallel.me==0)fprintf(stderr,"tau: %11.4e eII: %11.4e eta_old: %11.4e eta_new: %11.4e\n",tau, eedot[i],eta_old,eta_new);
}
}
}
@@ -1070,7 +1230,7 @@
ettnew*E->data.ref_viscosity);
#endif
// if(visited)
- // fprintf(stderr,"%11g %11g %11g %11g\n",ettnew,EEta[ (i-1)*vpts + jj ] ,ettby,eedot[i]);
+ // fprintf(stderr,"%11.4e %11.4e %11.4e %11.4e\n",ettnew,EEta[ (i-1)*vpts + jj ] ,ettby,eedot[i]);
EEta[ (i-1)*vpts + jj ] = ettnew;
}
} /* end regular plasticity */
@@ -1107,8 +1267,10 @@
/* determine composition of each of the nodes of the element */
for(kk = 1; kk <= ends; kk++){
CC[kk] = E->C[E->ien[i].node[kk]];
- if(CC[kk] < 0)CC[kk]=0.0;
- if(CC[kk] > 1)CC[kk]=1.0;
+ if(E->control.check_c_irange){
+ if(CC[kk] < 0)CC[kk]=0.0;
+ if(CC[kk] > 1)CC[kk]=1.0;
+ }
}
for(jj = 1; jj <= vpts; jj++)
{
@@ -1124,3 +1286,4 @@
}
visited++;
}
+
Modified: mc/3D/CitcomCU/trunk/src/anisotropic_viscosity.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/anisotropic_viscosity.h 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/anisotropic_viscosity.h 2010-10-11 04:21:32 UTC (rev 17256)
@@ -1,34 +1,48 @@
#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 ) ;
+#define CITCOM_ANIVISC_ORTHO_MODE 1
+#define CITCOM_ANIVISC_TI_MODE 2
+#define CITCOM_ANIVISC_ALIGN_WITH_VEL 0
+#define CITCOM_ANIVISC_ALIGN_WITH_ISA 1
+#define CITCOM_ANIVISC_MIXED_ALIGN 2
+
+void get_constitutive(double [6][6], int, int, double, double, int, struct All_variables *);
+void get_constitutive_ti_viscosity(double [6][6], double, double, double [3], int, double, double);
+void get_constitutive_orthotropic_viscosity(double [6][6], double, double [3], int, double, double);
+void get_constitutive_isotropic(double [6][6]);
+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 rotate_ti6x6_to_director(double [6][6], double [3]);
void get_citcom_spherical_rot(double, double, double [3][3]);
void get_orth_delta(double [3][3][3][3], double [3]);
+void align_director_with_ISA_for_element(struct All_variables *, int);
+unsigned char calc_isa_from_vgm(double [3][3], double [3], int, double [3], struct All_variables *, int);
+int is_pure_shear(double [3][3], double [3][3], double [3][3]);
void rot_4x4(double [3][3][3][3], double [3][3], double [3][3][3][3]);
+void zero_6x6(double [6][6]);
+void zero_4x4(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 print_6x6_mat(FILE *, 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 *);
+void isacalc(double [3][3], double *, double [3], struct All_variables *, int *);
+void f_times_ft(double [3][3], double [3][3]);
+void drex_eigen(double [3][3], double [3][3], int *);
+void malmul_scaled_id(double [3][3], double [3][3], double, double);
+
+void print_3x3_mat(FILE *, double [3][3]);
+
+void calc_exp_matrixt(double [3][3],double ,double [3][3],
+ struct All_variables *);
+
+void dgpadm_(int *,int *,double *,double *,int *,double *,int *,
+ int *,int *,int *,int *);
+
#define __CITCOM_READ_ANIVISC_HEADER__
#endif
Added: mc/3D/CitcomCU/trunk/src/expokit/dgpadm.f
===================================================================
--- mc/3D/CitcomCU/trunk/src/expokit/dgpadm.f (rev 0)
+++ mc/3D/CitcomCU/trunk/src/expokit/dgpadm.f 2010-10-11 04:21:32 UTC (rev 17256)
@@ -0,0 +1,170 @@
+*----------------------------------------------------------------------|
+ subroutine DGPADM( ideg,m,t,H,ldh,wsp,lwsp,ipiv,iexph,ns,iflag )
+
+ implicit none
+ integer ideg, m, ldh, lwsp, iexph, ns, iflag, ipiv(m)
+ double precision t, H(ldh,m), wsp(lwsp)
+
+*-----Purpose----------------------------------------------------------|
+*
+* Computes exp(t*H), the matrix exponential of a general matrix in
+* full, using the irreducible rational Pade approximation to the
+* exponential function exp(x) = r(x) = (+/-)( I + 2*(q(x)/p(x)) ),
+* combined with scaling-and-squaring.
+*
+*-----Arguments--------------------------------------------------------|
+*
+* ideg : (input) the degre of the diagonal Pade to be used.
+* a value of 6 is generally satisfactory.
+*
+* m : (input) order of H.
+*
+* H(ldh,m) : (input) argument matrix.
+*
+* t : (input) time-scale (can be < 0).
+*
+* wsp(lwsp) : (workspace/output) lwsp .ge. 4*m*m+ideg+1.
+*
+* ipiv(m) : (workspace)
+*
+*>>>> iexph : (output) number such that wsp(iexph) points to exp(tH)
+* i.e., exp(tH) is located at wsp(iexph ... iexph+m*m-1)
+* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+* NOTE: if the routine was called with wsp(iptr),
+* then exp(tH) will start at wsp(iptr+iexph-1).
+*
+* ns : (output) number of scaling-squaring used.
+*
+* iflag : (output) exit flag.
+* 0 - no problem
+* <0 - problem
+*
+*----------------------------------------------------------------------|
+* Roger B. Sidje (rbs at maths.uq.edu.au)
+* EXPOKIT: Software Package for Computing Matrix Exponentials.
+* ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
+*----------------------------------------------------------------------|
+*
+ integer mm,i,j,k,ih2,ip,iq,iused,ifree,iodd,icoef,iput,iget
+ double precision hnorm,scale,scale2,cp,cq
+
+ intrinsic INT,ABS,DBLE,LOG,MAX
+
+*--- check restrictions on input parameters ...
+ mm = m*m
+ iflag = 0
+ if ( ldh.lt.m ) iflag = -1
+ if ( lwsp.lt.4*mm+ideg+1 ) iflag = -2
+ if ( iflag.ne.0 ) stop 'bad sizes (in input of DGPADM)'
+*
+*--- initialise pointers ...
+*
+ icoef = 1
+ ih2 = icoef + (ideg+1)
+ ip = ih2 + mm
+ iq = ip + mm
+ ifree = iq + mm
+*
+*--- scaling: seek ns such that ||t*H/2^ns|| < 1/2;
+* and set scale = t/2^ns ...
+*
+ do i = 1,m
+ wsp(i) = 0.0d0
+ enddo
+ do j = 1,m
+ do i = 1,m
+ wsp(i) = wsp(i) + ABS( H(i,j) )
+ enddo
+ enddo
+ hnorm = 0.0d0
+ do i = 1,m
+ hnorm = MAX( hnorm,wsp(i) )
+ enddo
+ hnorm = ABS( t*hnorm )
+ if ( hnorm.eq.0.0d0 ) stop 'Error - null H in input of DGPADM.'
+ ns = MAX( 0,INT(LOG(hnorm)/LOG(2.0d0))+2 )
+ scale = t / DBLE(2**ns)
+ scale2 = scale*scale
+*
+*--- compute Pade coefficients ...
+*
+ i = ideg+1
+ j = 2*ideg+1
+ wsp(icoef) = 1.0d0
+ do k = 1,ideg
+ wsp(icoef+k) = (wsp(icoef+k-1)*DBLE( i-k ))/DBLE( k*(j-k) )
+ enddo
+*
+*--- H2 = scale2*H*H ...
+*
+ call DGEMM( 'n','n',m,m,m,scale2,H,ldh,H,ldh,0.0d0,wsp(ih2),m )
+*
+*--- initialize p (numerator) and q (denominator) ...
+*
+ cp = wsp(icoef+ideg-1)
+ cq = wsp(icoef+ideg)
+ do j = 1,m
+ do i = 1,m
+ wsp(ip + (j-1)*m + i-1) = 0.0d0
+ wsp(iq + (j-1)*m + i-1) = 0.0d0
+ enddo
+ wsp(ip + (j-1)*(m+1)) = cp
+ wsp(iq + (j-1)*(m+1)) = cq
+ enddo
+*
+*--- Apply Horner rule ...
+*
+ iodd = 1
+ k = ideg - 1
+ 100 continue
+ iused = iodd*iq + (1-iodd)*ip
+ call DGEMM( 'n','n',m,m,m, 1.0d0,wsp(iused),m,
+ . wsp(ih2),m, 0.0d0,wsp(ifree),m )
+ do j = 1,m
+ wsp(ifree+(j-1)*(m+1)) = wsp(ifree+(j-1)*(m+1))+wsp(icoef+k-1)
+ enddo
+ ip = (1-iodd)*ifree + iodd*ip
+ iq = iodd*ifree + (1-iodd)*iq
+ ifree = iused
+ iodd = 1-iodd
+ k = k-1
+ if ( k.gt.0 ) goto 100
+*
+*--- Obtain (+/-)(I + 2*(p\q)) ...
+*
+ if ( iodd .eq. 1 ) then
+ call DGEMM( 'n','n',m,m,m, scale,wsp(iq),m,
+ . H,ldh, 0.0d0,wsp(ifree),m )
+ iq = ifree
+ else
+ call DGEMM( 'n','n',m,m,m, scale,wsp(ip),m,
+ . H,ldh, 0.0d0,wsp(ifree),m )
+ ip = ifree
+ endif
+ call DAXPY( mm, -1.0d0,wsp(ip),1, wsp(iq),1 )
+ call DGESV( m,m, wsp(iq),m, ipiv, wsp(ip),m, iflag )
+ if ( iflag.ne.0 ) stop 'Problem in DGESV (within DGPADM)'
+ call DSCAL( mm, 2.0d0, wsp(ip), 1 )
+ do j = 1,m
+ wsp(ip+(j-1)*(m+1)) = wsp(ip+(j-1)*(m+1)) + 1.0d0
+ enddo
+ iput = ip
+ if ( ns.eq.0 .and. iodd.eq.1 ) then
+ call DSCAL( mm, -1.0d0, wsp(ip), 1 )
+ goto 200
+ endif
+*
+*-- squaring : exp(t*H) = (exp(t*H))^(2^ns) ...
+*
+ iodd = 1
+ do k = 1,ns
+ iget = iodd*ip + (1-iodd)*iq
+ iput = (1-iodd)*ip + iodd*iq
+ call DGEMM( 'n','n',m,m,m, 1.0d0,wsp(iget),m, wsp(iget),m,
+ . 0.0d0,wsp(iput),m )
+ iodd = 1-iodd
+ enddo
+ 200 continue
+ iexph = iput
+ END
+*----------------------------------------------------------------------|
Modified: mc/3D/CitcomCU/trunk/src/global_defs.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/global_defs.h 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/global_defs.h 2010-10-11 04:21:32 UTC (rev 17256)
@@ -43,6 +43,7 @@
#ifdef USE_GGRD
#include "hc.h"
+
#endif
#if defined(__osf__)
@@ -757,10 +758,10 @@
float Ra_410, clapeyron410, transT410, width410;
#ifdef USE_GGRD
struct ggrd_master ggrd;
- int slab_slice;
- float slab_theta_bound;
+ int ggrd_slab_slice;
+ float ggrd_slab_theta_bound[4];
+ struct ggrd_gt ggrd_ss_grd[4];
-
#endif
@@ -985,6 +986,7 @@
float *VIn1[MAX_LEVELS],*EVIn1[MAX_LEVELS];
float *VIn2[MAX_LEVELS],*EVIn2[MAX_LEVELS];
float *VIn3[MAX_LEVELS],*EVIn3[MAX_LEVELS];
+ unsigned char *avmode[MAX_LEVELS];
#endif
@@ -1041,7 +1043,20 @@
/* To regenerate function prototypes:
* 1. Comment out the #include line below, and save this file
- * 2. Run command: cproto -q -p -f 3 *.c > prototypes.h
+ * 2. Run command: cproto -q -p -f 2 *.c > prototypes.h
* 3. Uncomment the #include line below, and save this file
+
+ could also use
+ cproto -DCITCOM_ALLOW_ANISOTROPIC_VISC -DUSE_GGRD -q -p -I. -f2 *.c > prototypes.h
+
*/
+void twiddle_thumbs(struct All_variables *, int ); /* somehow, cproto
+ does not pick up
+ these functions?! */
+void convection_initial_temperature_ggrd(struct All_variables *);
+void set_2dc_defaults(struct All_variables *);
+void remove_horiz_ave(struct All_variables *, float *, float *, int );
+void output_velo_related_gzdir(struct All_variables *, int);
+void output_velo_related(struct All_variables *, int);
+
#include "prototypes.h"
Modified: mc/3D/CitcomCU/trunk/src/prototypes.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/prototypes.h 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/prototypes.h 2010-10-11 04:21:32 UTC (rev 17256)
@@ -11,6 +11,33 @@
void std_timestep(struct All_variables *);
void process_heating(struct All_variables *);
/* Anisotropic_viscosity.c */
+void get_constitutive(double [6][6], int, int, double, double, int, struct All_variables *);
+void get_constitutive_ti_viscosity(double [6][6], double, double, double [3], int, double, double);
+void get_constitutive_orthotropic_viscosity(double [6][6], double, double [3], int, double, double);
+void get_constitutive_isotropic(double [6][6]);
+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 rotate_ti6x6_to_director(double [6][6], double [3]);
+void get_citcom_spherical_rot(double, double, double [3][3]);
+void get_orth_delta(double [3][3][3][3], double [3]);
+void align_director_with_ISA_for_element(struct All_variables *, int);
+unsigned char calc_isa_from_vgm(double [3][3], double [3], int, double [3], struct All_variables *, int);
+int is_pure_shear(double [3][3], double [3][3], double [3][3]);
+void rot_4x4(double [3][3][3][3], double [3][3], double [3][3][3][3]);
+void zero_6x6(double [6][6]);
+void zero_4x4(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 print_6x6_mat(FILE *, double [6][6]);
+void c4fromc6(double [3][3][3][3], double [6][6]);
+void c6fromc4(double [6][6], double [3][3][3][3]);
+void drex_isacalc(double [3][3], double *, double [3], struct All_variables *, int *);
+void f_times_ft(double [3][3], double [3][3]);
+void drex_eigen(double [3][3], double [3][3], int *);
+void malmul_scaled_id(double [3][3], double [3][3], double, double);
/* Boundary_conditions.c */
void velocity_boundary_conditions(struct All_variables *);
void temperature_boundary_conditions(struct All_variables *);
@@ -76,6 +103,9 @@
/* Element_calculations.c */
void assemble_forces(struct All_variables *, int);
void get_elt_k(struct All_variables *, int, double [24 * 24], int, int);
+void get_ba(struct Shape_function *, struct Shape_function_dx *, struct CC *, struct CCX *, double [4][9], int, int, int, int, double [9][4][9][7]);
+void get_ba_p(struct Shape_function *, struct Shape_function_dx *, struct CC *, struct CCX *, double [4][9], int, int, int, int, double [9][4][9][7]);
+void get_vgm_p(double [4][9], struct Shape_function *, struct Shape_function_dx *, struct CC *, struct CCX *, double [4][9], int, int, int, int, double [3][3], double [3]);
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);
@@ -112,7 +142,7 @@
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);
+int solve_del2_u(struct All_variables *, double *, double *, double, 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);
@@ -134,10 +164,13 @@
void set_3ds_defaults(struct All_variables *);
void set_3dc_defaults(struct All_variables *);
/* Ggrd_handling.c */
+void ggrd_solve_eigen3x3(double [3][3], double [3], double [3][3], struct All_variables *);
+void ggrd_read_anivisc_from_file(struct All_variables *);
/* Global_operations.c */
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_div_norm2(struct All_variables *, double *);
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);
@@ -209,7 +242,18 @@
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 xyz2rtp(float, float, float, float *);
void myerror(char *, struct All_variables *);
+void get_9vec_from_3x3(double *, double [3][3]);
+void get_3x3_from_9vec(double [3][3], double *);
+void get_3x3_from_6vec(double [3][3], double *);
+void mat_mult_vec_3x3(double [3][3], double [3], double [3]);
+void normalize_vec3(float *, float *, float *);
+void normalize_vec3d(double *, double *, double *);
+void cross_product(double [3], double [3], double [3]);
+void matmul_3x3(double [3][3], double [3][3], double [3][3]);
+void assign_to_3x3(double [3][3], double);
+void remove_trace_3x3(double [3][3]);
/* Parallel_related.c */
void parallel_process_termination(void);
void parallel_domain_decomp1(struct All_variables *);
@@ -243,11 +287,14 @@
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 *, float *, float *, float *, float *);
/* Process_buoyancy.c */
+void process_temp_field(struct All_variables *, int);
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 *, int);
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 *);
@@ -255,9 +302,11 @@
/* Profiling.c */
float CPU_time(void);
/* Shape_functions.c */
+void construct_shape_functions(struct All_variables *);
double lpoly(int, double);
double lpolydash(int, double);
/* Size_does_matter.c */
+void twiddle_thumbs(struct All_variables *, int);
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);
@@ -266,6 +315,7 @@
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 *);
void cg_allocate_vars(struct All_variables *);
void assemble_forces_iterative(struct All_variables *);
/* Solver_multigrid.c */
@@ -282,6 +332,7 @@
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 *);
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 *);
@@ -289,7 +340,10 @@
void compute_sphereh_table(struct All_variables *);
/* Stokes_flow_Incomp.c */
void solve_constrained_flow_iterative(struct All_variables *);
+float solve_Ahat_p_fhat_new(struct All_variables *, double *, double *, double *, double, int *);
+void initial_vel_residual(struct All_variables *, double *, double *, double *, double *, double);
float solve_Ahat_p_fhat(struct All_variables *, double *, double *, double *, double, int *);
+void generate_log_message(int, double, double, double, double, struct All_variables *);
void v_from_vector(struct All_variables *, float **, double *);
/* Topo_gravity.c */
void get_CBF_topo(struct All_variables *, float *, float *);
@@ -304,15 +358,13 @@
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);
+double second_invariant_from_3x3(double [3][3]);
+void calc_vgm_matrix(struct All_variables *, double *, double *);
+void calc_strain_from_vgm(double [3][3], double [3][3]);
+void calc_strain_from_vgm9(double *, double [3][3]);
+void calc_rot_from_vgm(double [3][3], double [3][3]);
+void calc_strain_rate_matrix(struct All_variables *, double *);
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]);
-
+int in_slab_slice(float , int , struct All_variables *);
Modified: mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h 2010-10-10 23:20:29 UTC (rev 17255)
+++ mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h 2010-10-11 04:21:32 UTC (rev 17256)
@@ -39,6 +39,8 @@
viscosity fields, those determined from prior input, those
related to temperature/pressure/stress/anything else. */
+#define CITCOM_CU_VISC_MAXLAYER 40
+
struct VISC_OPT
{
void (*update_viscosity) ();
@@ -54,9 +56,10 @@
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 */
+ int anisotropic_init; /* 0: isotropic, 1: random, 2: from file 3: align with ISA */
char anisotropic_init_dir[1000];
int anivisc_layer; /* layer to assign anisotropic viscosity to for mode 2 */
+ double ani_vis2_factor; /* for mode 3, anisotropy scale factor*/
#endif
@@ -77,7 +80,7 @@
float zlith;
float zcomp;
- float zbase_layer[40]; /* new */
+ float zbase_layer[CITCOM_CU_VISC_MAXLAYER]; /* new */
int FREEZE;
float freeze_thresh;
@@ -92,14 +95,14 @@
float sdepv_misfit;
float sdepv_iter_damp;
int sdepv_normalize;
- float sdepv_expt[40];
- float sdepv_trns[40];
+ float sdepv_expt[CITCOM_CU_VISC_MAXLAYER];
+ float sdepv_trns[CITCOM_CU_VISC_MAXLAYER];
int TDEPV;
int TDEPV_AVE;
- float N0[40];
- float E[40], T0[40];
- float T[40], Z[40];
+ float N0[CITCOM_CU_VISC_MAXLAYER];
+ float E[CITCOM_CU_VISC_MAXLAYER], T0[CITCOM_CU_VISC_MAXLAYER];
+ float T[CITCOM_CU_VISC_MAXLAYER], Z[CITCOM_CU_VISC_MAXLAYER];
/* byerlee :
@@ -117,8 +120,8 @@
lithostatic pressure. Caution: 0 might produce crap.
*/
int BDEPV;
- float abyerlee[40],bbyerlee[40],
- lbyerlee[40];
+ float abyerlee[CITCOM_CU_VISC_MAXLAYER],bbyerlee[CITCOM_CU_VISC_MAXLAYER],
+ lbyerlee[CITCOM_CU_VISC_MAXLAYER];
float plasticity_viscosity_offset;
@@ -136,29 +139,29 @@
int CDEPV; /* composition dependent viscosity */
- float pre_comp[40]; /* prefactors */
+ float pre_comp[CITCOM_CU_VISC_MAXLAYER]; /* prefactors */
int weak_blobs;
- float weak_blobx[40];
- float weak_bloby[40];
- float weak_blobz[40];
- float weak_blobwidth[40];
- float weak_blobmag[40];
+ float weak_blobx[CITCOM_CU_VISC_MAXLAYER];
+ float weak_bloby[CITCOM_CU_VISC_MAXLAYER];
+ float weak_blobz[CITCOM_CU_VISC_MAXLAYER];
+ float weak_blobwidth[CITCOM_CU_VISC_MAXLAYER];
+ float weak_blobmag[CITCOM_CU_VISC_MAXLAYER];
int weak_zones;
- float weak_zonex1[40];
- float weak_zoney1[40];
- float weak_zonez1[40];
- float weak_zonex2[40];
- float weak_zoney2[40];
- float weak_zonez2[40];
+ float weak_zonex1[CITCOM_CU_VISC_MAXLAYER];
+ float weak_zoney1[CITCOM_CU_VISC_MAXLAYER];
+ float weak_zonez1[CITCOM_CU_VISC_MAXLAYER];
+ float weak_zonex2[CITCOM_CU_VISC_MAXLAYER];
+ float weak_zoney2[CITCOM_CU_VISC_MAXLAYER];
+ float weak_zonez2[CITCOM_CU_VISC_MAXLAYER];
- float weak_zonewidth[40];
- float weak_zonemag[40];
+ float weak_zonewidth[CITCOM_CU_VISC_MAXLAYER];
+ float weak_zonemag[CITCOM_CU_VISC_MAXLAYER];
int guess;
char old_file[100];
@@ -168,8 +171,7 @@
char VISC_OPT[20];
int layers; /* number of layers with properties .... */
- float layer_depth[40];
- float layer_visc[40];
+ float layer_depth[CITCOM_CU_VISC_MAXLAYER];
int SLABLVZ; /* slab structure imposed on top of 3 layer structure */
int slvzd1, slvzd2, slvzd3; /* layer thicknesses (nodes) */
@@ -187,10 +189,10 @@
/* MODULE BASED VISCOSITY VARIATIONS */
int RESDEPV;
- float RESeta0[40];
+ float RESeta0[CITCOM_CU_VISC_MAXLAYER];
int CHEMDEPV;
- float CH0[40];
- float CHEMeta0[40];
+ float CH0[CITCOM_CU_VISC_MAXLAYER];
+ float CHEMeta0[CITCOM_CU_VISC_MAXLAYER];
} viscosity;
More information about the CIG-COMMITS
mailing list