[cig-commits] r8124 - in mc/3D/CitcomS/trunk: bin examples/Full lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Oct 17 11:35:04 PDT 2007
Author: tan2
Date: 2007-10-17 11:35:03 -0700 (Wed, 17 Oct 2007)
New Revision: 8124
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:
Roll back r8111, reading velo files for init. temperature stays as tic_method=-1. Other changes in r8111 will be put back later.
Modified: mc/3D/CitcomS/trunk/bin/Citcom.c
===================================================================
--- mc/3D/CitcomS/trunk/bin/Citcom.c 2007-10-17 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/bin/Citcom.c 2007-10-17 18:35:03 UTC (rev 8124)
@@ -68,7 +68,7 @@
float dot();
float cpu_time_on_vp_it;
- int cpu_total_seconds,k, *temp,out_cycles;
+ int cpu_total_seconds,k, *temp;
double CPU_time0(),time,initial_time,start_time,avaimem();
struct All_variables *E;
@@ -88,7 +88,6 @@
start_time = time = CPU_time0();
read_instructions(E, argv[1]);
-
initial_setup(E);
cpu_time_on_vp_it = CPU_time0();
@@ -101,16 +100,16 @@
}
/* This if-block is replaced by CitcomS.Solver.launch()*/
- if ((E->control.restart == 1) || E->control.post_p) {
+ if (E->control.restart || E->control.post_p) {
read_checkpoint(E);
if (E->control.post_p) {
post_processing(E);
parallel_process_termination();
}
- }else {
- /* regular init, or restart from T only */
- initial_conditions(E);
+ }
+ else {
+ initial_conditions(E);
if(E->control.pseudo_free_surf) {
if(E->mesh.topvbc == 2)
@@ -122,23 +121,12 @@
general_stokes_solver(E);
}
- 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);
- }
+ (E->problem_output)(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)
@@ -182,28 +170,20 @@
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 ((out_cycles % E->output.write_q_files)==0)
+ if ((E->monitor.solution_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, out_cycles);
+ output_time(E, E->monitor.solution_cycles);
- if ((out_cycles % E->control.checkpoint_frequency)==0)
- output_checkpoint(E);
+ if ((E->monitor.solution_cycles % E->control.checkpoint_frequency)==0) {
+ output_checkpoint(E);
+ }
if(E->control.mat_control==1)
read_mat_from_file(E);
@@ -233,13 +213,13 @@
if (E->parallel.me == 0) {
- fprintf(stderr,"cycles=%d\n",out_cycles);
+ fprintf(stderr,"cycles=%d\n",E->monitor.solution_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 ? (1):(0)))));
+ cpu_time_on_vp_it/((float)(E->monitor.solution_cycles-E->control.restart)));
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 ? (1):(0)))));
+ cpu_time_on_vp_it/((float)(E->monitor.solution_cycles-E->control.restart)));
}
output_finalize(E);
Modified: mc/3D/CitcomS/trunk/examples/Full/input.sample
===================================================================
--- mc/3D/CitcomS/trunk/examples/Full/input.sample 2007-10-17 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/examples/Full/input.sample 2007-10-17 18:35:03 UTC (rev 8124)
@@ -90,7 +90,7 @@
# miscellaneous information
stokes_flow_only=0
-inputdiffusivity=1.0
+inputdiffusicity=0.00100
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-17 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c 2007-10-17 18:35:03 UTC (rev 8124)
@@ -109,7 +109,6 @@
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)];
@@ -127,7 +126,6 @@
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)];
@@ -151,8 +149,6 @@
/* 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 */
@@ -173,21 +169,15 @@
omega[i1] += lomega[i*3+i1] * vtmp;
vw += vtmp;
}
- 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;
+ for(i1=0;i1 < 3;i1++)
+ omega[i1] /= vw;
+
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;
}
@@ -226,7 +216,7 @@
amp = 0.0;
break;
case 1: /* add this velocity */
- if((fabs(theta) > 1e-5) &&(fabs(theta-M_PI) > 1e-5)){
+ if((fabs(theta) > 1e-6) &&(fabs(theta-M_PI)>1e-6)){
coslat=sin(theta);
coslon=cos(phi);
sinlat=cos(theta);
@@ -242,8 +232,7 @@
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;
@@ -252,7 +241,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;
@@ -273,7 +262,6 @@
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);
@@ -289,7 +277,7 @@
//
// subtract a net rotation component from a velocity
-// field given as v_theta (velt) and v_phi (velp)
+// field given as v_theta (velt) and v_phi (velp) on
//
void sub_netr(float r,float theta,float phi,float *velt,float *velp, double *omega)
@@ -324,8 +312,8 @@
vphi = vx * px + vy * py;
/* remove */
- *velt = *velt - vtheta;
- *velp = *velp - vphi;
+ *velt -= vtheta;
+ *velp -= vphi;
}
Modified: mc/3D/CitcomS/trunk/lib/Drive_solvers.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-10-17 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-10-17 18:35:03 UTC (rev 8124)
@@ -68,7 +68,7 @@
float vmag;
- double Udot_mag, dUdot_mag,omega[3];
+ double Udot_mag, dUdot_mag;
int m,count,i,j,k;
double *oldU[NCS], *delta_U[NCS];
@@ -104,8 +104,7 @@
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++) {
@@ -124,13 +123,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*/
@@ -143,9 +142,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-17 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-10-17 18:35:03 UTC (rev 8124)
@@ -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-17 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-10-17 18:35:03 UTC (rev 8124)
@@ -32,11 +32,6 @@
#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
@@ -791,16 +786,12 @@
return;
}
-/*
-remove rigid rotation every time step
- */
-
void remove_rigid_rot(struct All_variables *E)
{
-
- double myatan();
+ void velo_from_element_d();
+ double myatan();
double wx, wy, wz, v_theta, v_phi;
double vx[9], vy[9], vz[9];
double r, t, f;
@@ -814,8 +805,7 @@
const int ppts = PPOINTS3D;
const int vpts = VPOINTS3D;
const int sphere_key = 1;
-
- float VV[4][9];
+ double VV[4][9];
double rot, fr, tr;
/* Note: no need to weight in rho(r) here. */
@@ -832,7 +822,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];
@@ -845,12 +835,11 @@
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];
@@ -873,7 +862,6 @@
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-17 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2007-10-17 18:35:03 UTC (rev 8124)
@@ -81,8 +81,7 @@
*/
switch(E->convection.tic_method){
- case -1: /* restart from file, no other options
- needed */
+ case -1: /* restart from file */
break;
case 0:
case 3:
@@ -190,8 +189,7 @@
if (E->control.lith_age)
lith_age_construct_tic(E);
- else if (E->control.restart == 2) {
- /* read restart info for T from file */
+ else if (E->convection.tic_method == -1) {
#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-17 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-10-17 18:35:03 UTC (rev 8124)
@@ -84,7 +84,6 @@
void initial_mesh_solver_setup(struct All_variables *E)
{
-
E->monitor.cpu_time_at_last_cycle =
E->monitor.cpu_time_at_start = CPU_time0();
@@ -184,6 +183,7 @@
================================================== */
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,13 +526,11 @@
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);
@@ -540,7 +538,6 @@
}
-
/* ===================================
Functions which set up details
common to all problems follow ...
@@ -1082,7 +1079,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
@@ -1103,7 +1100,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
@@ -1253,13 +1250,14 @@
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();
- }
+ 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();
+ }
}
}
@@ -1357,7 +1355,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-17 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2007-10-17 18:35:03 UTC (rev 8124)
@@ -74,12 +74,7 @@
*/
input_int("gzdir_vtkio",&(E->output.gzdir.vtk_io),"0",m);
- input_boolean("gzdir_rnr",&(E->output.gzdir.rnr),"off",m); /* remove
- net
- rotation
- on
- output
- step? */
+ input_boolean("gzdir_rnr",&(E->output.gzdir.rnr),"off",m); /* remove net rotation? */
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-17 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-10-17 18:35:03 UTC (rev 8124)
@@ -127,12 +127,12 @@
/**********************************************************************/
-void gzdir_output(struct All_variables *E, int out_cycles)
+void gzdir_output(struct All_variables *E, int cycles)
{
char output_dir[255];
- if (out_cycles == 0 ){
+ if (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,out_cycles);
+ snprintf(output_dir,255,"%s/%d",E->control.data_dir,cycles);
mkdatadir(output_dir);
/* output */
- gzdir_output_velo_temp(E, out_cycles); /* don't move this around,
+ gzdir_output_velo_temp(E, cycles); /* don't move this around,
else new VTK output won't
work */
- gzdir_output_visc(E, out_cycles);
+ gzdir_output_visc(E, cycles);
- gzdir_output_surf_botm(E, out_cycles);
+ gzdir_output_surf_botm(E, cycles);
/* optiotnal output below */
/* compute and output geoid (in spherical harmonics coeff) */
if (E->output.geoid)
- gzdir_output_geoid(E, out_cycles);
+ gzdir_output_geoid(E, cycles);
if (E->output.stress)
- gzdir_output_stress(E, out_cycles);
+ gzdir_output_stress(E, cycles);
if (E->output.pressure)
- gzdir_output_pressure(E, out_cycles);
+ gzdir_output_pressure(E, cycles);
if (E->output.horiz_avg)
- gzdir_output_horiz_avg(E, out_cycles);
+ gzdir_output_horiz_avg(E, cycles);
if(E->control.tracer){
if(E->output.tracer ||
- (out_cycles == E->advection.max_timesteps))
- gzdir_output_tracer(E, out_cycles);
+ (cycles == E->advection.max_timesteps))
+ gzdir_output_tracer(E, cycles);
}
if (E->output.comp_nd && E->composition.on)
- gzdir_output_comp_nd(E, out_cycles);
+ gzdir_output_comp_nd(E, cycles);
if (E->output.comp_el && E->composition.on)
- gzdir_output_comp_el(E, out_cycles);
+ gzdir_output_comp_el(E, cycles);
if(E->output.heating && E->control.disptn_number != 0)
- gzdir_output_heating(E, out_cycles);
+ gzdir_output_heating(E, cycles);
return;
}
@@ -470,9 +470,6 @@
}
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",
@@ -550,13 +547,9 @@
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]; /* 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]
- */
+ vcorr[0] = E->sphere.cap[j].V[1][i];
+ vcorr[1] = E->sphere.cap[j].V[2][i];
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;
@@ -1171,11 +1164,7 @@
}
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++) {
@@ -1184,10 +1173,6 @@
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;
}
}
@@ -1212,7 +1197,6 @@
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 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-10-17 18:35:03 UTC (rev 8124)
@@ -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,sq_vdotv;
+ double alpha, delta, r0dotz0, r1dotz1;
double residual, v_res;
double global_vdot(), global_pdot();
@@ -206,18 +206,15 @@
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 %6.2f s v=%.3e div/v=%.3e "
+ fprintf(E->fp,"AhatP (%03d) after %g seconds with div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
+ count, CPU_time0()-time0, E->monitor.incompressibility,
0.0, 0.0, E->monitor.solution_cycles);
fflush(E->fp);
- fprintf(stderr,"AhatP (%03d) after %6.2f s v=%.3e div/v=%.3e "
+ fprintf(stderr,"AhatP (%03d) after %g seconds with div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
+ count, CPU_time0()-time0, E->monitor.incompressibility,
0.0, 0.0, E->monitor.solution_cycles);
}
@@ -307,16 +304,14 @@
count++;
- sq_vdotv = sqrt(E->monitor.vdotv);
-
if(E->control.print_convergence && E->parallel.me==0) {
- fprintf(E->fp, "AhatP (%03d) after %6.2f s v=%.3e div/v=%.3e "
+ fprintf(E->fp, "AhatP (%03d) after %g seconds with div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
+ count, CPU_time0()-time0, E->monitor.incompressibility,
dvelocity, dpressure, E->monitor.solution_cycles);
- fprintf(stderr, "AhatP (%03d) after %6.2f s v=%.3e div/v=%.3e "
+ fprintf(stderr, "AhatP (%03d) after %g seconds with div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, sq_vdotv, E->monitor.incompressibility,
+ count, CPU_time0()-time0, E->monitor.incompressibility,
dvelocity, dpressure, E->monitor.solution_cycles);
}
@@ -374,7 +369,7 @@
int m, j, count, lev;
int valid;
- double alpha, beta, omega,sq_vdotv;
+ double alpha, beta, omega;
double r0dotrt, r1dotrt;
double residual, dpressure, dvelocity;
@@ -433,18 +428,15 @@
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 s, v=%.3e div/v=%.3e "
+ fprintf(E->fp,"AhatP (%03d) after %g seconds with div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
+ count, CPU_time0()-time0, E->monitor.incompressibility,
0.0, 0.0, E->monitor.solution_cycles);
fflush(E->fp);
- fprintf(stderr,"AhatP (%03d) after %g s, v=%.3e div/v=%.3e "
+ fprintf(stderr,"AhatP (%03d) after %g seconds with div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
+ count, CPU_time0()-time0, E->monitor.incompressibility,
0.0, 0.0, E->monitor.solution_cycles);
}
@@ -573,17 +565,15 @@
count++;
-
- sq_vdotv = sqrt(E->monitor.vdotv);
if(E->control.print_convergence && E->parallel.me==0) {
- fprintf(E->fp, "AhatP (%03d) after %g s, v=%.3e, div/v=%.3e "
+ fprintf(E->fp, "AhatP (%03d) after %g seconds with div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, sq_vdotv ,E->monitor.incompressibility,
+ count, CPU_time0()-time0, E->monitor.incompressibility,
dvelocity, dpressure, E->monitor.solution_cycles);
- fprintf(stderr, "AhatP (%03d) after %g s, v=%.3e div/v=%.3e "
+ fprintf(stderr, "AhatP (%03d) after %g seconds with div/v=%.3e "
"dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
+ count, CPU_time0()-time0, 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 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-10-17 18:35:03 UTC (rev 8124)
@@ -42,7 +42,6 @@
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)
@@ -84,7 +83,7 @@
E->viscosity.sdepv_misfit = 1.0;
input_boolean("SDEPV",&(E->viscosity.SDEPV),"off",m);
if (E->viscosity.SDEPV) {
- E->viscosity.sdepv_visited = 0;
+ E->viscosity.pdepv_visited = 0;
input_float_vector("sdepv_expt",E->viscosity.num_mat,(E->viscosity.sdepv_expt),m);
}
@@ -105,12 +104,9 @@
input_boolean("CDEPV",&(E->viscosity.CDEPV),"off",m);
- 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->viscosity.CDEPV){ /* compositional viscosity */
+ if(!E->control.tracer)
+ myerror(E,"error: CDEPV requires tracers, but tracer is off");
if(E->trace.nflavors > 10)
myerror(E,"error: too many flavors for CDEPV");
/* read in flavor factors */
@@ -534,20 +530,13 @@
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];
@@ -613,10 +602,10 @@
}else{
for(e=1;e<=nel;e++) /* initialize with unity if no velocities around */
- eedot[e] = 1.0;
+ eedot[e] = 1.0e-5;
if(m == E->sphere.caps_per_proc)
E->viscosity.pdepv_visited = 1;
- if((E->parallel.me == 0)&&(E->control.verbose)){
+ if(E->parallel.me == 0){
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 00:28:12 UTC (rev 8123)
+++ mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2007-10-17 18:35:03 UTC (rev 8124)
@@ -81,7 +81,7 @@
int SDEPV;
float sdepv_misfit;
- int sdepv_normalize,sdepv_visited;
+ int sdepv_normalize;
float sdepv_expt[40];
float sdepv_trns[40];
More information about the cig-commits
mailing list