[cig-commits] r8131 - mc/3D/CitcomS/trunk/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Oct 17 13:55:17 PDT 2007
Author: tan2
Date: 2007-10-17 13:55:17 -0700 (Wed, 17 Oct 2007)
New Revision: 8131
Modified:
mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c
mc/3D/CitcomS/trunk/lib/Drive_solvers.c
mc/3D/CitcomS/trunk/lib/Output_gzdir.c
mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
Log:
Put the rest of r8111 back
Modified: mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c 2007-10-17 20:33:21 UTC (rev 8130)
+++ mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c 2007-10-17 20:55:17 UTC (rev 8131)
@@ -109,6 +109,7 @@
x[0] += E->x[m][1][lnode[d]] * vtmp; /* coords */
x[1] += E->x[m][2][lnode[d]] * vtmp;
x[2] += E->x[m][3][lnode[d]] * vtmp;
+
v[0] += E->sphere.cap[m].V[1][lnode[d]] * vtmp; /* theta */
v[1] += E->sphere.cap[m].V[2][lnode[d]] * vtmp; /* phi */
vw += dGamma.vpt[GMVGAMMA(0,nint)];
@@ -126,6 +127,7 @@
x[0] += E->x[m][1][lnode[d]] * vtmp; /* coords */
x[1] += E->x[m][2][lnode[d]] * vtmp;
x[2] += E->x[m][3][lnode[d]] * vtmp;
+ /* */
v[0] += E->sphere.cap[m].V[1][lnode[d]] * vtmp;
v[1] += E->sphere.cap[m].V[2][lnode[d]] * vtmp;
vw += dGamma.vpt[GMVGAMMA(1,nint)];
@@ -149,6 +151,8 @@
/* depth range */
rr = E->sx[1][3][E->ien[1][elz].node[5]] - E->sx[1][3][E->ien[1][1].node[1]];
+ if(rr < 1e-7)
+ myerror(E,"rr error in net r determine");
vw = 0.0;
for (i=0;i < elz;i++) { /* regular 0..n-1 loop */
/* solve layer NR */
@@ -169,15 +173,21 @@
omega[i1] += lomega[i*3+i1] * vtmp;
vw += vtmp;
}
- for(i1=0;i1 < 3;i1++)
- omega[i1] /= vw;
-
+ if(fabs(vw) > 1e-8) /* when would it be zero? */
+ for(i1=0;i1 < 3;i1++)
+ omega[i1] /= vw;
+ else
+ for(i1=0;i1 < 3;i1++)
+ omega[i1] = 0.0;
free ((void *) acoef);
free ((void *) coef);
free ((void *) lomega);
oamp = sqrt(omega[0]*omega[0] + omega[1]*omega[1] + omega[2]*omega[2]);
+ if(E->parallel.me == 0)
+ fprintf(stderr,"determined net rotation of | %.4e %.4e %.4e | = %.4e\n",
+ omega[0],omega[1],omega[2],oamp);
return oamp;
}
@@ -216,7 +226,7 @@
amp = 0.0;
break;
case 1: /* add this velocity */
- if((fabs(theta) > 1e-6) &&(fabs(theta-M_PI)>1e-6)){
+ if((fabs(theta) > 1e-5) &&(fabs(theta-M_PI) > 1e-5)){
coslat=sin(theta);
coslon=cos(phi);
sinlat=cos(theta);
@@ -232,7 +242,8 @@
b = -rz*coslon;
c = rz*rzu*coslon+rx*coslat;
d = -rz*sinlon;
- e = -ry*rzu*coslon+rx*rzu*sinlon;
+ e = -ry*rzu*coslon+rx*rzu*sinlon;
+
f = ry*sinlon+rx*coslon;
c9[0] += a*a+b*b;
@@ -241,7 +252,7 @@
c9[3] += c*c+d*d;
c9[4] += c*e+d*f;
c9[5] += e*e+f*f;
-
+
c9[6] += a*velt+b*velp;
c9[7] += c*velt+d*velp;
c9[8] += e*velt+f*velp;
@@ -262,6 +273,7 @@
omega[0] = c9[6];
omega[1] = c9[7];
omega[2] = c9[8];
+
/* solve solution*/
hc_ludcmp_3x3(coef,ind);
hc_lubksb_3x3(coef,ind,omega);
@@ -277,7 +289,7 @@
//
// subtract a net rotation component from a velocity
-// field given as v_theta (velt) and v_phi (velp) on
+// field given as v_theta (velt) and v_phi (velp)
//
void sub_netr(float r,float theta,float phi,float *velt,float *velp, double *omega)
@@ -312,8 +324,8 @@
vphi = vx * px + vy * py;
/* remove */
- *velt -= vtheta;
- *velp -= vphi;
+ *velt = *velt - vtheta;
+ *velp = *velp - vphi;
}
Modified: mc/3D/CitcomS/trunk/lib/Drive_solvers.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-10-17 20:33:21 UTC (rev 8130)
+++ mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-10-17 20:55:17 UTC (rev 8131)
@@ -68,7 +68,7 @@
float vmag;
- double Udot_mag, dUdot_mag;
+ double Udot_mag, dUdot_mag,omega[3];
int m,count,i,j,k;
double *oldU[NCS], *delta_U[NCS];
@@ -104,7 +104,8 @@
Udot_mag=dUdot_mag=0.0;
count=1;
- while (1) {
+ while (1) {
+
for (m=1;m<=E->sphere.caps_per_proc;m++)
for (i=0;i<neq;i++) {
@@ -123,13 +124,13 @@
dUdot_mag,Udot_mag,count);
fflush(E->fp);
}
- if ((count>50) || (dUdot_mag<E->viscosity.sdepv_misfit))
+ if ((count>50) || (dUdot_mag < E->viscosity.sdepv_misfit))
break;
-
+
get_system_viscosity(E,1,E->EVI[E->mesh.levmax],E->VI[E->mesh.levmax]);
construct_stiffness_B_matrix(E);
solve_constrained_flow_iterative(E);
-
+
count++;
} /*end while*/
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-10-17 20:33:21 UTC (rev 8130)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-10-17 20:55:17 UTC (rev 8131)
@@ -127,12 +127,12 @@
/**********************************************************************/
-void gzdir_output(struct All_variables *E, int cycles)
+void gzdir_output(struct All_variables *E, int out_cycles)
{
char output_dir[255];
- if (cycles == 0) {
+ if (out_cycles == 0 ){
/* initial I/O */
-
+
gzdir_output_coord(E);
/*gzdir_output_mat(E);*/
}
@@ -145,48 +145,48 @@
*/
/* make a directory */
- snprintf(output_dir,255,"%s/%d",E->control.data_dir,cycles);
+ snprintf(output_dir,255,"%s/%d",E->control.data_dir,out_cycles);
mkdatadir(output_dir);
/* output */
- gzdir_output_velo_temp(E, cycles); /* don't move this around,
+ gzdir_output_velo_temp(E, out_cycles); /* don't move this around,
else new VTK output won't
work */
- gzdir_output_visc(E, cycles);
+ gzdir_output_visc(E, out_cycles);
- gzdir_output_surf_botm(E, cycles);
+ gzdir_output_surf_botm(E, out_cycles);
/* optiotnal output below */
/* compute and output geoid (in spherical harmonics coeff) */
if (E->output.geoid)
- gzdir_output_geoid(E, cycles);
+ gzdir_output_geoid(E, out_cycles);
if (E->output.stress)
- gzdir_output_stress(E, cycles);
+ gzdir_output_stress(E, out_cycles);
if (E->output.pressure)
- gzdir_output_pressure(E, cycles);
+ gzdir_output_pressure(E, out_cycles);
if (E->output.horiz_avg)
- gzdir_output_horiz_avg(E, cycles);
+ gzdir_output_horiz_avg(E, out_cycles);
if(E->control.tracer){
if(E->output.tracer ||
- (cycles == E->advection.max_timesteps))
- gzdir_output_tracer(E, cycles);
+ (out_cycles == E->advection.max_timesteps))
+ gzdir_output_tracer(E, out_cycles);
}
if (E->output.comp_nd && E->composition.on)
- gzdir_output_comp_nd(E, cycles);
+ gzdir_output_comp_nd(E, out_cycles);
if (E->output.comp_el && E->composition.on)
- gzdir_output_comp_el(E, cycles);
+ gzdir_output_comp_el(E, out_cycles);
if(E->output.heating && E->control.disptn_number != 0)
- gzdir_output_heating(E, cycles);
+ gzdir_output_heating(E, out_cycles);
return;
}
@@ -470,6 +470,9 @@
}
if(E->output.gzdir.rnr){ /* remove the whole model net rotation */
+ if((E->control.remove_rigid_rotation)&&
+ (E->parallel.me == 0)) /* that's not too terrible but wastes time */
+ fprintf(stderr,"WARNING: both gzdir.rnr and remove_rigid_rotation are switched on!\n");
oamp = determine_model_net_rotation(E,omega);
if(E->parallel.me == 0)
fprintf(stderr,"gzdir_output_velo_temp: removing net rotation: |%8.3e, %8.3e, %8.3e| = %8.3e\n",
@@ -547,9 +550,13 @@
if(E->output.gzdir.rnr){
/* remove NR */
for(i=1;i<=E->lmesh.nno;i++,k += 9) {
- vcorr[0] = E->sphere.cap[j].V[1][i];
- vcorr[1] = E->sphere.cap[j].V[2][i];
+ vcorr[0] = E->sphere.cap[j].V[1][i]; /* vtheta */
+ vcorr[1] = E->sphere.cap[j].V[2][i]; /* vphi */
+ /* remove the velocity that corresponds to a net rotation of omega[0..2] at location
+ r,t,p from the t,p velocities in vcorr[0..1]
+ */
sub_netr(E->sx[j][3][i],E->sx[j][1][i],E->sx[j][2][i],(vcorr+0),(vcorr+1),omega);
+
convert_pvec_to_cvec(E->sphere.cap[j].V[3][i],vcorr[0],vcorr[1],
(E->output.gzdir.vtk_base+k),cvec);
if(be_write_float_to_file(cvec,3,fp1)!=3)BE_WERROR;
@@ -1164,7 +1171,11 @@
}
if(fscanf(fp,"%i %i %f",&ll,&mm,&restart_elapsed_time) != 3)
myerror(E,"restart vtkl read error 0");
-
+ if(mm != E->lmesh.nno){
+ fprintf(stderr,"%i %i\n",mm, E->lmesh.nno);
+ myerror(E,"lmesh.nno mismatch in restart files");
+ }
+
switch(E->output.gzdir.vtk_io) {
case 1: /* VTK */
for(m=1;m <= E->sphere.caps_per_proc;m++) {
@@ -1173,6 +1184,10 @@
for(i=1;i<=E->lmesh.nno;i++){
if(fscanf(fp,"%f",&g) != 1)
myerror(E,"restart vtkl read error 2");
+ if(!finite(g)){
+ fprintf(stderr,"WARNING: found a NaN in input temperatures\n");
+ g=0.0;
+ }
E->T[m][i] = g;
}
}
@@ -1197,6 +1212,7 @@
gzip_file(output_file);
temperatures_conform_bcs(E);
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-10-17 20:33:21 UTC (rev 8130)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-10-17 20:55:17 UTC (rev 8131)
@@ -137,7 +137,7 @@
double *r1[NCS], *r2[NCS], *z1[NCS], *s1[NCS], *s2[NCS], *F[NCS];
double *shuffle[NCS];
- double alpha, delta, r0dotz0, r1dotz1;
+ double alpha, delta, r0dotz0, r1dotz1,sq_vdotv;
double residual, v_res;
double global_vdot(), global_pdot();
@@ -206,15 +206,18 @@
residual = incompressibility_residual(E, V, r1);
+
+ sq_vdotv = sqrt(E->monitor.vdotv);
+
if (E->control.print_convergence && E->parallel.me==0) {
- fprintf(E->fp,"AhatP (%03d) after %g seconds with div/v=%.3e "
+ fprintf(E->fp,"AhatP (%03d) after %6.2f s v=%.3e div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, E->monitor.incompressibility,
+ count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
0.0, 0.0, E->monitor.solution_cycles);
fflush(E->fp);
- fprintf(stderr,"AhatP (%03d) after %g seconds with div/v=%.3e "
+ fprintf(stderr,"AhatP (%03d) after %6.2f s v=%.3e div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, E->monitor.incompressibility,
+ count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
0.0, 0.0, E->monitor.solution_cycles);
}
@@ -304,14 +307,16 @@
count++;
+ sq_vdotv = sqrt(E->monitor.vdotv);
+
if(E->control.print_convergence && E->parallel.me==0) {
- fprintf(E->fp, "AhatP (%03d) after %g seconds with div/v=%.3e "
+ fprintf(E->fp, "AhatP (%03d) after %6.2f s v=%.3e div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, E->monitor.incompressibility,
+ count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
dvelocity, dpressure, E->monitor.solution_cycles);
- fprintf(stderr, "AhatP (%03d) after %g seconds with div/v=%.3e "
+ fprintf(stderr, "AhatP (%03d) after %6.2f s v=%.3e div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, E->monitor.incompressibility,
+ count, CPU_time0()-time0, sq_vdotv, E->monitor.incompressibility,
dvelocity, dpressure, E->monitor.solution_cycles);
}
@@ -369,7 +374,7 @@
int m, j, count, lev;
int valid;
- double alpha, beta, omega;
+ double alpha, beta, omega,sq_vdotv;
double r0dotrt, r1dotrt;
double residual, dpressure, dvelocity;
@@ -428,15 +433,18 @@
rt[m][j] = r1[m][j];
+ sq_vdotv = sqrt(E->monitor.vdotv);
+
+
if (E->control.print_convergence && E->parallel.me==0) {
- fprintf(E->fp,"AhatP (%03d) after %g seconds with div/v=%.3e "
+ fprintf(E->fp,"AhatP (%03d) after %g s, v=%.3e div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, E->monitor.incompressibility,
+ count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
0.0, 0.0, E->monitor.solution_cycles);
fflush(E->fp);
- fprintf(stderr,"AhatP (%03d) after %g seconds with div/v=%.3e "
+ fprintf(stderr,"AhatP (%03d) after %g s, v=%.3e div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, E->monitor.incompressibility,
+ count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
0.0, 0.0, E->monitor.solution_cycles);
}
@@ -565,15 +573,17 @@
count++;
+
+ sq_vdotv = sqrt(E->monitor.vdotv);
if(E->control.print_convergence && E->parallel.me==0) {
- fprintf(E->fp, "AhatP (%03d) after %g seconds with div/v=%.3e "
+ fprintf(E->fp, "AhatP (%03d) after %g s, v=%.3e, div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, E->monitor.incompressibility,
+ count, CPU_time0()-time0, sq_vdotv ,E->monitor.incompressibility,
dvelocity, dpressure, E->monitor.solution_cycles);
- fprintf(stderr, "AhatP (%03d) after %g seconds with div/v=%.3e "
+ fprintf(stderr, "AhatP (%03d) after %g s, v=%.3e div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, E->monitor.incompressibility,
+ count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
dvelocity, dpressure, E->monitor.solution_cycles);
}
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-10-17 20:33:21 UTC (rev 8130)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-10-17 20:55:17 UTC (rev 8131)
@@ -42,6 +42,7 @@
static void apply_low_visc_wedge_channel(struct All_variables *E, float **evisc);
static void low_viscosity_channel_factor(struct All_variables *E, float *F);
static void low_viscosity_wedge_factor(struct All_variables *E, float *F);
+void parallel_process_termination();
void viscosity_system_input(struct All_variables *E)
@@ -83,7 +84,7 @@
E->viscosity.sdepv_misfit = 1.0;
input_boolean("SDEPV",&(E->viscosity.SDEPV),"off",m);
if (E->viscosity.SDEPV) {
- E->viscosity.pdepv_visited = 0;
+ E->viscosity.sdepv_visited = 0;
input_float_vector("sdepv_expt",E->viscosity.num_mat,(E->viscosity.sdepv_expt),m);
}
@@ -104,9 +105,12 @@
input_boolean("CDEPV",&(E->viscosity.CDEPV),"off",m);
- if(E->viscosity.CDEPV){ /* compositional viscosity */
- if(!E->control.tracer)
- myerror(E,"error: CDEPV requires tracers, but tracer is off");
+ if(E->viscosity.CDEPV){
+ /* compositional viscosity */
+ if(E->control.tracer < 1){
+ fprintf(stderr,"error: CDEPV requires tracers, but tracer is off\n");
+ parallel_process_termination();
+ }
if(E->trace.nflavors > 10)
myerror(E,"error: too many flavors for CDEPV");
/* read in flavor factors */
@@ -530,13 +534,20 @@
two = 2.0;
for(m=1;m<=E->sphere.caps_per_proc;m++) {
-
+ if(E->viscosity.sdepv_visited){
+
/* get second invariant for all elements */
strain_rate_2_inv(E,m,eedot,1);
+ }else{
+ for(e=1;e<=nel;e++) /* initialize with unity if no velocities around */
+ eedot[e] = 1.0;
+ E->viscosity.sdepv_visited = 1;
+ }
/* eedot cannot be too small, or the viscosity will go to inf */
- for(e=1;e<=nel;e++)
- eedot[e] = max(eedot[e], 1.0e-16);
+ for(e=1;e<=nel;e++){
+ eedot[e] = max(eedot[e], 1.0e-16);
+ }
for(e=1;e<=nel;e++) {
exponent1= one/E->viscosity.sdepv_expt[E->mat[m][e]-1];
@@ -602,10 +613,10 @@
}else{
for(e=1;e<=nel;e++) /* initialize with unity if no velocities around */
- eedot[e] = 1.0e-5;
+ eedot[e] = 1.0;
if(m == E->sphere.caps_per_proc)
E->viscosity.pdepv_visited = 1;
- if(E->parallel.me == 0){
+ if((E->parallel.me == 0)&&(E->control.verbose)){
for(e=0;e < E->viscosity.num_mat;e++)
fprintf(stderr,"num mat: %i a: %g b: %g y: %g\n",
e,E->viscosity.pdepv_a[e],E->viscosity.pdepv_b[e],E->viscosity.pdepv_y[e]);
Modified: mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2007-10-17 20:33:21 UTC (rev 8130)
+++ mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2007-10-17 20:55:17 UTC (rev 8131)
@@ -81,7 +81,7 @@
int SDEPV;
float sdepv_misfit;
- int sdepv_normalize;
+ int sdepv_normalize,sdepv_visited;
float sdepv_expt[40];
float sdepv_trns[40];
More information about the cig-commits
mailing list