[cig-commits] r8111 - in mc/3D/CitcomS/trunk: bin examples/Full lib
becker at geodynamics.org
becker at geodynamics.org
Sat Oct 13 12:07:54 PDT 2007
Author: becker
Date: 2007-10-13 12:07:53 -0700 (Sat, 13 Oct 2007)
New Revision: 8111
Modified:
mc/3D/CitcomS/trunk/bin/Citcom.c
mc/3D/CitcomS/trunk/examples/Full/input.sample
mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c
mc/3D/CitcomS/trunk/lib/Drive_solvers.c
mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
mc/3D/CitcomS/trunk/lib/Global_operations.c
mc/3D/CitcomS/trunk/lib/Initial_temperature.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Output.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:
Fixed rigid rotation code bug.
Added restart=2 option which will not use checkpoints, but tic_method=-1 style temperature input
(formally realized via tic_method=-1, but this should be cleaner)
fixed input diffusivity parameter in full example file.
Modified: mc/3D/CitcomS/trunk/bin/Citcom.c
===================================================================
--- mc/3D/CitcomS/trunk/bin/Citcom.c 2007-10-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/bin/Citcom.c 2007-10-13 19:07:53 UTC (rev 8111)
@@ -68,7 +68,7 @@
float dot();
float cpu_time_on_vp_it;
- int cpu_total_seconds,k, *temp;
+ int cpu_total_seconds,k, *temp,out_cycles;
double CPU_time0(),time,initial_time,start_time,avaimem();
struct All_variables *E;
@@ -88,6 +88,7 @@
start_time = time = CPU_time0();
read_instructions(E, argv[1]);
+
initial_setup(E);
cpu_time_on_vp_it = CPU_time0();
@@ -100,16 +101,16 @@
}
/* This if-block is replaced by CitcomS.Solver.launch()*/
- if (E->control.restart || E->control.post_p) {
+ if ((E->control.restart == 1) || E->control.post_p) {
read_checkpoint(E);
if (E->control.post_p) {
post_processing(E);
parallel_process_termination();
}
- }
- else {
- initial_conditions(E);
+ }else {
+ /* regular init, or restart from T only */
+ initial_conditions(E);
if(E->control.pseudo_free_surf) {
if(E->mesh.topvbc == 2)
@@ -121,12 +122,23 @@
general_stokes_solver(E);
}
- (E->problem_output)(E, E->monitor.solution_cycles);
+ if(strcmp(E->output.format, "ascii-gz") == 0){ /* gzdir output likes continuous cycle numbering
+ in output files didn't want to change elsewhere
+ since sometimes cycles == 0
+ is used to test for inits TWB
+ */
+ if(!E->control.restart) { /* no restart */
+ /* first step */
+ (E->problem_output)(E, E->monitor.solution_cycles);
+ output_time(E, E->monitor.solution_cycles);
+ }
+ }else{
+ /* regular output */
+ (E->problem_output)(E, E->monitor.solution_cycles);
+ /* information about simulation time and wall clock time */
+ output_time(E, E->monitor.solution_cycles);
+ }
- /* information about simulation time and wall clock time */
- output_time(E, E->monitor.solution_cycles);
-
-
if (E->control.stokes) {
if(E->control.tracer==1)
@@ -170,20 +182,28 @@
tracer_advection(E);
general_stokes_solver(E);
+
+ /* the output cycle counter is cumulative for VTK IO TWB */
+ if(strcmp(E->output.format, "ascii-gz") == 0) /* vtk output */
+ if(E->control.restart) /* restart */
+ out_cycles = E->monitor.solution_cycles + E->monitor.solution_cycles_init;
+ else
+ out_cycles = E->monitor.solution_cycles ;
+ else /* other output */
+ out_cycles = E->monitor.solution_cycles ;
+
+
if(E->output.write_q_files)
- if ((E->monitor.solution_cycles % E->output.write_q_files)==0)
+ if ((out_cycles % E->output.write_q_files)==0)
heat_flux(E);
+ if ((out_cycles % E->control.record_every)==0)
+ (E->problem_output)(E, out_cycles);
- if ((E->monitor.solution_cycles % E->control.record_every)==0) {
- (E->problem_output)(E, E->monitor.solution_cycles);
- }
-
/* information about simulation time and wall clock time */
- output_time(E, E->monitor.solution_cycles);
+ output_time(E, out_cycles);
- if ((E->monitor.solution_cycles % E->control.checkpoint_frequency)==0) {
- output_checkpoint(E);
- }
+ if ((out_cycles % E->control.checkpoint_frequency)==0)
+ output_checkpoint(E);
if(E->control.mat_control==1)
read_mat_from_file(E);
@@ -213,13 +233,13 @@
if (E->parallel.me == 0) {
- fprintf(stderr,"cycles=%d\n",E->monitor.solution_cycles);
+ fprintf(stderr,"cycles=%d\n",out_cycles);
cpu_time_on_vp_it=CPU_time0()-cpu_time_on_vp_it;
fprintf(stderr,"Average cpu time taken for velocity step = %f\n",
- cpu_time_on_vp_it/((float)(E->monitor.solution_cycles-E->control.restart)));
+ cpu_time_on_vp_it/((float)(E->monitor.solution_cycles-(E->control.restart ? (1):(0)))));
fprintf(E->fp,"Initialization overhead = %f\n",initial_time);
fprintf(E->fp,"Average cpu time taken for velocity step = %f\n",
- cpu_time_on_vp_it/((float)(E->monitor.solution_cycles-E->control.restart)));
+ cpu_time_on_vp_it/((float)(E->monitor.solution_cycles-(E->control.restart ? (1):(0)))));
}
output_finalize(E);
Modified: mc/3D/CitcomS/trunk/examples/Full/input.sample
===================================================================
--- mc/3D/CitcomS/trunk/examples/Full/input.sample 2007-10-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/examples/Full/input.sample 2007-10-13 19:07:53 UTC (rev 8111)
@@ -90,7 +90,7 @@
# miscellaneous information
stokes_flow_only=0
-inputdiffusicity=0.00100
+inputdiffusivity=1.0
rayleigh=5.0e06
Q0=0
Modified: mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c 2007-10-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c 2007-10-13 19:07:53 UTC (rev 8111)
@@ -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-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-10-13 19:07:53 UTC (rev 8111)
@@ -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*/
@@ -142,9 +143,9 @@
} /*end if SDEPV or PDEPV */
/* remove the rigid rotation component from the velocity solution */
- if(E->sphere.caps == 12 && E->control.remove_rigid_rotation)
- remove_rigid_rot(E);
-
+ if((E->sphere.caps == 12) && E->control.remove_rigid_rotation){
+ remove_rigid_rot(E);
+ }
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-10-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-10-13 19:07:53 UTC (rev 8111)
@@ -368,7 +368,7 @@
}
}
break;
- case -1:
+ case -1:
case 1:
case 2:
break;
Modified: mc/3D/CitcomS/trunk/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-10-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-10-13 19:07:53 UTC (rev 8111)
@@ -32,6 +32,11 @@
#include "element_definitions.h"
#include "global_defs.h"
+
+void velo_from_element(struct All_variables *,float [4][9],int ,int ,int);
+
+
+
/* ===============================================
strips horizontal average from nodal field X.
Assumes orthogonal mesh, otherwise, horizontals
@@ -786,12 +791,16 @@
return;
}
+/*
+remove rigid rotation every time step
+ */
+
void remove_rigid_rot(struct All_variables *E)
{
- void velo_from_element_d();
- double myatan();
+
+ double myatan();
double wx, wy, wz, v_theta, v_phi;
double vx[9], vy[9], vz[9];
double r, t, f;
@@ -805,7 +814,8 @@
const int ppts = PPOINTS3D;
const int vpts = VPOINTS3D;
const int sphere_key = 1;
- double VV[4][9];
+
+ float VV[4][9];
double rot, fr, tr;
/* Note: no need to weight in rho(r) here. */
@@ -822,7 +832,7 @@
for (e=1;e<=E->lmesh.nel;e++) {
- t = E->eco[m][e].centre[1];
+ t = E->eco[m][e].centre[1];
f = E->eco[m][e].centre[2];
r = E->eco[m][e].centre[3];
@@ -835,11 +845,12 @@
for (j=1;j<=ppts;j++) {
for (i=1;i<=ends;i++) {
+
vx[j] += VV[1][i]*E->N.ppt[GNPINDEX(i,j)];
vy[j] += VV[2][i]*E->N.ppt[GNPINDEX(i,j)];
}
}
-
+
wx = -r*vy[1];
wy = r*vx[1];
@@ -862,6 +873,7 @@
if (E->parallel.me==0) {
fprintf(E->fp,"Rigid rotation: rot=%e tr=%e fr=%e\n",rot,tr*180/M_PI,fr*180/M_PI);
+ fprintf(stderr,"Rigid rotation: rot=%e tr=%e fr=%e\n",rot,tr*180/M_PI,fr*180/M_PI);
}
Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2007-10-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2007-10-13 19:07:53 UTC (rev 8111)
@@ -81,7 +81,8 @@
*/
switch(E->convection.tic_method){
- case -1: /* restart from file */
+ case -1: /* restart from file, no other options
+ needed */
break;
case 0:
case 3:
@@ -189,7 +190,8 @@
if (E->control.lith_age)
lith_age_construct_tic(E);
- else if (E->convection.tic_method == -1) {
+ else if (E->control.restart == 2) {
+ /* read restart info for T from file */
#ifdef USE_GZDIR
if(strcmp(E->output.format, "ascii-gz") == 0)
restart_tic_from_gzdir_file(E);
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-10-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-10-13 19:07:53 UTC (rev 8111)
@@ -84,6 +84,7 @@
void initial_mesh_solver_setup(struct All_variables *E)
{
+
E->monitor.cpu_time_at_last_cycle =
E->monitor.cpu_time_at_start = CPU_time0();
@@ -183,7 +184,6 @@
================================================== */
setup_parser(E,filename);
-
global_default_values(E);
read_initial_settings(E);
@@ -385,7 +385,7 @@
input_boolean("see_convergence",&(E->control.print_convergence),"off",m);
input_int("stokes_flow_only",&(E->control.stokes),"0",m);
-
+ /* 1: restart from checkpoint file 2: restart from tic style file */
input_int("restart",&(E->control.restart),"0",m);
input_int("post_p",&(E->control.post_p),"0",m);
input_int("solution_cycles_init",&(E->monitor.solution_cycles_init),"0",m);
@@ -526,11 +526,13 @@
lith_age_input(E);
tic_input(E);
+
tracer_input(E);
viscosity_input(E); /* moved the viscosity input behind
the tracer input */
+
(E->problem_settings)(E);
@@ -538,6 +540,7 @@
}
+
/* ===================================
Functions which set up details
common to all problems follow ...
@@ -1079,7 +1082,7 @@
else
sprintf(logfile,"%s.log", E->control.data_file);
- if (E->control.restart || E->control.post_p)
+ if (E->control.restart || E->control.post_p )
/* append the log file if restart */
E->fp = output_open(logfile, "a");
else
@@ -1100,7 +1103,7 @@
else
sprintf(timeoutput,"%s.time", E->control.data_file);
- if (E->control.restart || E->control.post_p)
+ if (E->control.restart || E->control.post_p )
/* append the time file if restart */
E->fptime = output_open(timeoutput, "a");
else
@@ -1250,14 +1253,13 @@
parallel_process_termination();
}
- if (E->control.restart || E->control.post_p ||
- E->convection.tic_method == -1 ||
- (E->control.tracer && E->trace.ic_method == 2)) {
- found = strchr(E->control.data_prefix_old, '/');
- if (found) {
- fprintf(stderr, "error in input parameter: datafile_old='%s' contains '/'\n", E->control.data_file);
- parallel_process_termination();
- }
+ if ((E->control.restart==1) || E->control.post_p ||
+ (E->control.tracer && (E->trace.ic_method == 2))) {
+ found = strchr(E->control.data_prefix_old, '/');
+ if (found) {
+ fprintf(stderr, "error in input parameter: datafile_old='%s' contains '/'\n", E->control.data_file);
+ parallel_process_termination();
+ }
}
}
@@ -1355,7 +1357,7 @@
E->control.data_prefix);
if (E->control.restart || E->control.post_p ||
- E->convection.tic_method == -1 ||
+ (E->convection.tic_method == -1) ||
(E->control.tracer && E->trace.ic_method == 2)) {
expand_datadir(E, E->control.data_dir_old);
snprintf(E->control.old_P_file, 200, "%s/%s", E->control.data_dir_old,
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2007-10-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2007-10-13 19:07:53 UTC (rev 8111)
@@ -74,7 +74,12 @@
*/
input_int("gzdir_vtkio",&(E->output.gzdir.vtk_io),"0",m);
- input_boolean("gzdir_rnr",&(E->output.gzdir.rnr),"off",m); /* remove net rotation? */
+ input_boolean("gzdir_rnr",&(E->output.gzdir.rnr),"off",m); /* remove
+ net
+ rotation
+ on
+ output
+ step? */
E->output.gzdir.vtk_base_init = 0;
E->output.gzdir.vtk_base_save = 1; /* should we save the basis vectors? (memory!) */
//fprintf(stderr,"gzdir: vtkio: %i save basis vectors: %i\n",
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-10-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-10-13 19:07:53 UTC (rev 8111)
@@ -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-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-10-13 19:07:53 UTC (rev 8111)
@@ -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-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-10-13 19:07:53 UTC (rev 8111)
@@ -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-12 20:24:17 UTC (rev 8110)
+++ mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2007-10-13 19:07:53 UTC (rev 8111)
@@ -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