[cig-commits] r7790 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Wed Aug 8 18:27:33 PDT 2007
Author: becker
Date: 2007-08-08 18:27:32 -0700 (Wed, 08 Aug 2007)
New Revision: 7790
Modified:
mc/3D/CitcomS/trunk/lib/Initial_temperature.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Output_gzdir.c
mc/3D/CitcomS/trunk/lib/Problem_related.c
mc/3D/CitcomS/trunk/lib/Tracer_setup.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
Log:
- added a compositional viscosity function, CDEPV, based on two tracer
flavors
- for this to work, I had to move viscosity_input() *behind*
tic_input() and tracer_input() in Instructions.c
- built in restart facilities for temperature and tracers when using
gzdir I/O
Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2007-08-09 00:36:01 UTC (rev 7789)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2007-08-09 01:27:32 UTC (rev 7790)
@@ -40,6 +40,9 @@
void debug_tic(struct All_variables *);
void restart_tic_from_file(struct All_variables *);
+#ifdef USE_GZDIR
+void restart_tic_from_gzdir_file(struct All_variables *);
+#endif
void tic_input(struct All_variables *E)
@@ -166,7 +169,12 @@
lith_age_restart_tic(E);
else
*/
- restart_tic_from_file(E);
+#ifdef USE_GZDIR
+ if(strcmp(E->output.format, "gzdir") == 0)
+ restart_tic_from_gzdir_file(E);
+ else
+#endif
+ restart_tic_from_file(E);
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-09 00:36:01 UTC (rev 7789)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-09 01:27:32 UTC (rev 7790)
@@ -398,10 +398,13 @@
h5input_params(E);
phase_change_input(E);
lith_age_input(E);
- viscosity_input(E);
+
tic_input(E);
tracer_input(E);
+ viscosity_input(E); /* moved the viscosity input behind
+ the tracer input */
+
(E->problem_settings)(E);
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-08-09 00:36:01 UTC (rev 7789)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-08-09 01:27:32 UTC (rev 7790)
@@ -65,11 +65,15 @@
void gzdir_output_pressure(struct All_variables *, int);
+void restart_tic_from_gzdir_file(struct All_variables *);
+
void calc_cbase_at_tp(float , float , float *);
void rtp2xyz(float , float , float, float *);
void convert_pvec_to_cvec(float ,float , float , float *,float *);
void *safe_malloc (size_t );
+int open_file_zipped(char *, FILE **,struct All_variables *);
+void gzip_file(char *);
extern void parallel_process_termination();
@@ -127,13 +131,15 @@
if (E->output.horiz_avg == 1)
gzdir_output_horiz_avg(E, cycles);
- if(E->output.tracer == 1 && E->control.tracer == 1)
+ if(E->control.tracer){
+ if(E->output.tracer || (cycles == E->advection.max_timesteps))
gzdir_output_tracer(E, cycles);
+ }
- if (E->output.comp_nd == 1 && E->composition.on)
+ if (E->output.comp_nd && E->composition.on)
gzdir_output_comp_nd(E, cycles);
- if (E->output.comp_el == 1 && E->composition.on)
+ if (E->output.comp_el && E->composition.on)
gzdir_output_comp_el(E, cycles);
return;
@@ -544,7 +550,8 @@
E->parallel.me, cycles);
fp1 = gzdir_output_open(output_file,"w");
- gzprintf(fp1,"%d %d %.5e\n",cycles,E->lmesh.nno,E->monitor.elapsed_time);
+ gzprintf(fp1,"%d %d %.5e\n",cycles,E->lmesh.nno,
+ E->monitor.elapsed_time);
for(j=1;j<=E->sphere.caps_per_proc;j++) {
gzprintf(fp1,"%3d %7d\n",j,E->lmesh.nno);
@@ -565,9 +572,9 @@
char output_file[255];
gzFile *fp1;
- snprintf(output_file,255,"%s/%d/tracer.%d.%d.gz", E->control.data_dir,
- cycles,
- E->parallel.me, cycles);
+ snprintf(output_file,255,"%s/%d/tracer.%d.%d.gz",
+ E->control.data_dir,cycles,
+ E->parallel.me, cycles);
fp1 = gzdir_output_open(output_file,"w");
ncolumns = 3 + E->trace.number_of_extra_quantities;
@@ -617,7 +624,8 @@
E->composition.bulk_composition);
for(i=1;i<=E->lmesh.nno;i++) {
- gzprintf(fp1,"%.6e\n",E->composition.comp_node[j][i]);
+ gzprintf(fp1,"%.6e\n",
+ E->composition.comp_node[j][i]);
}
}
@@ -653,7 +661,127 @@
return;
}
+/*
+restart facility for zipped/VTK style , will init temperature
+*/
+void restart_tic_from_gzdir_file(struct All_variables *E)
+{
+ int ii, ll, mm,rezip;
+ float restart_elapsed_time;
+ int i, m;
+ char output_file[255], input_s[1000];
+ FILE *fp;
+ float v1, v2, v3, g;
+
+ ii = E->monitor.solution_cycles_init;
+
+ if(E->output.gzdir_vtkio) { /* VTK I/O */
+ snprintf(output_file,255,"%s/%d/t.%d.%d",
+ E->control.data_dir_old,
+ ii,E->parallel.me,ii);
+ }else{ /* old output */
+ snprintf(output_file,255,"%s/%d/velo.%d.%d",
+ E->control.data_dir_old,ii,
+ E->parallel.me,ii);
+ }
+ /* open file */
+ rezip = open_file_zipped(output_file,&fp,E);
+ if (E->parallel.me==0)
+ fprintf(E->fp,"restart_tic_from_gzdir_file: using %s for restarted temperature\n",
+ output_file);
+
+ if(fscanf(fp,"%i %i %f",&ll,&mm,&restart_elapsed_time) != 3)
+ myerror(E,"restart vtkl read error 0");
+
+ if(E->output.gzdir_vtkio) { /* VTK */
+ for(m=1;m <= E->sphere.caps_per_proc;m++) {
+ if(fscanf(fp,"%i %i",&ll,&mm) != 2)
+ myerror(E,"restart vtkl read error 1");
+ for(i=1;i<=E->lmesh.nno;i++)
+ if(fscanf(fp,"%f",&(E->T[m][i]))!=1)
+ myerror(E,"restart vtkl read error 2");
+ }
+ }else{ /* old style velo */
+ for(m=1;m <= E->sphere.caps_per_proc;m++) {
+ fscanf(fp,"%i %i",&ll,&mm);
+ for(i=1;i<=E->lmesh.nno;i++) {
+ fscanf(fp,"%f %f %f %f",&v1,&v2,&v3,&g);
+ /* E->sphere.cap[m].V[1][i] = v1;
+ E->sphere.cap[m].V[1][i] = v2;
+ E->sphere.cap[m].V[1][i] = v3; */
+ //E->T[m][i] = max(0.0,min(g,1.0));
+ E->T[m][i] = g;
+ }
+ }
+ }
+ fclose (fp);
+ if(rezip) /* rezip */
+ gzip_file(output_file);
+
+ temperatures_conform_bcs(E);
+ return;
+}
+
+
+/*
+
+tries to open 'name'. if name exists, out will be pointer to file and
+return 0. if name doesn't exist, will check for name.gz. if this
+exists, will unzip and open, and return 1
+
+the idea is to preserve the initial file state
+
+*/
+int open_file_zipped(char *name, FILE **in,
+ struct All_variables *E)
+{
+ char mstring[1000];
+ *in = fopen(name,"r");
+ if (*in == NULL) {
+ /*
+ unzipped file not found
+ */
+ snprintf(mstring,1000,"%s.gz",name);
+ *in= fopen(mstring,"r");
+ if(*in != NULL){
+ /*
+ zipped version was found
+ */
+ fclose(*in);
+ snprintf(mstring,1000,"gunzip -f %s.gz",name); /* brutal */
+ system(mstring); /* unzip */
+ /* open unzipped file for read */
+ *in = fopen(name,"r");
+ if(*in == NULL)
+ myerror("open_file_zipped: unzipping error",E);
+ return 1;
+ }else{
+ /*
+ no file, either zipped or unzipped
+ */
+ snprintf(mstring,1000,"no files %s and %s.gz were found, exiting",
+ name,name);
+ myerror(E,mstring);
+ return 0;
+ }
+ }else{
+ /*
+ file was found unzipped
+ */
+ return 0;
+ }
+}
+
+/* compress a file using the sytem command */
+void gzip_file(char *output_file)
+{
+ char command_string[300];
+ snprintf(command_string,300,"gzip -f %s",output_file); /* brutal */
+ system(command_string);
+}
+
+
#endif /* gzdir switch */
Modified: mc/3D/CitcomS/trunk/lib/Problem_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Problem_related.c 2007-08-09 00:36:01 UTC (rev 7789)
+++ mc/3D/CitcomS/trunk/lib/Problem_related.c 2007-08-09 01:27:32 UTC (rev 7790)
@@ -30,6 +30,10 @@
#include "element_definitions.h"
#include "global_defs.h"
+#ifdef USE_GZDIR
+int open_file_zipped(char *, FILE **,struct All_variables *);
+void gzip_file(char *);
+#endif
/*=======================================================================
read velocity vectors at the top surface from files
=========================================================================*/
@@ -62,23 +66,43 @@
struct All_variables *E;
{
FILE *fp;
- int ll, mm;
+ int ll, mm,rezip;
char output_file[255],input_s[1000];
E->monitor.elapsed_time = 0.0;
if ((E->control.restart || E->control.post_p)) {
+
+#ifdef USE_GZDIR /* gzdir output */
+ if(strcmp(E->output.format, "gzdir") == 0){
+ if(E->output.gzdir_vtkio)
+ sprintf(output_file, "%s/%d/t.%d.%d",
+ E->control.data_dir_old,E->monitor.solution_cycles_init,E->parallel.me,E->monitor.solution_cycles_init);
+ else
+ sprintf(output_file, "%s/%d/velo.%d.%d",
+ E->control.data_dir_old,E->monitor.solution_cycles_init,E->parallel.me,E->monitor.solution_cycles_init);
+ }else{
sprintf(output_file, "%s.velo.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
- fp=fopen(output_file,"r");
- if (fp == NULL) {
- fprintf(E->fp,"(Problem_related #8) Cannot open %s\n",output_file);
- exit(8);
- }
- fgets(input_s,1000,fp);
- sscanf(input_s,"%d %d %f",&ll,&mm,&E->monitor.elapsed_time);
- fclose(fp);
- } /* end control.restart */
+ }
+ rezip = open_file_zipped(output_file,&fp,E);
+#else /* all others */
+ sprintf(output_file, "%s.velo.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
+ fp=fopen(output_file,"r");
+#endif
- return;
+ if (fp == NULL) {
+ fprintf(E->fp,"(Problem_related #8) Cannot open %s\n",output_file);
+ exit(8);
+ }
+ fgets(input_s,1000,fp);
+ sscanf(input_s,"%d %d %f",&ll,&mm,&E->monitor.elapsed_time);
+ fclose(fp);
+#ifdef USE_GZDIR
+ if(rezip)
+ gzip_file(output_file);
+#endif
+ } /* end control.restart */
+
+ return;
}
/*=======================================================================
Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-08-09 00:36:01 UTC (rev 7789)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-08-09 01:27:32 UTC (rev 7790)
@@ -48,6 +48,10 @@
#ifdef USE_GGRD
void ggrd_init_tracer_flavors(struct All_variables *);
#endif
+#ifdef USE_GZDIR
+int open_file_zipped(char *, FILE **,struct All_variables *);
+void gzip_file(char *);
+#endif
int icheck_that_processor_shell(struct All_variables *E,
int j, int nprocessor, double rad);
@@ -925,7 +929,7 @@
char output_file[200];
char input_s[1000];
- int i,j,kk;
+ int i,j,kk,rezip;
int idum1,ncolumns;
int numtracers;
@@ -944,13 +948,30 @@
parallel_process_termination();
}
+
+
+ /* deal with different output formats */
+#ifdef USE_GZDIR
+ if(strcmp(E->output.format, "gzdir") == 0){
+ sprintf(output_file,"%s/%d/tracer.%d.%d",
+ E->control.data_dir_old,E->monitor.solution_cycles_init,E->parallel.me,E->monitor.solution_cycles_init);
+ rezip = open_file_zipped(output_file,&fp1,E);
+ }else{
+ sprintf(output_file,"%s.tracer.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
+ if ( (fp1=fopen(output_file,"r"))==NULL) {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-file not found %s\n",output_file);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+#else
sprintf(output_file,"%s.tracer.%d.%d",E->control.old_P_file,E->parallel.me,E->monitor.solution_cycles_init);
-
if ( (fp1=fopen(output_file,"r"))==NULL) {
fprintf(E->trace.fpt,"ERROR(restart tracers)-file not found %s\n",output_file);
fflush(E->trace.fpt);
exit(10);
}
+#endif
fprintf(stderr,"Restarting Tracers from %s\n",output_file);
fflush(stderr);
@@ -1012,8 +1033,12 @@
}
fclose(fp1);
+#ifdef USE_GZDIR
+ if(strcmp(E->output.format, "gzdir") == 0)
+ if(rezip) /* rezip */
+ gzip_file(output_file);
+#endif
-
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-09 00:36:01 UTC (rev 7789)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-09 01:27:32 UTC (rev 7790)
@@ -59,7 +59,10 @@
}
+ for(i=0;i<10;i++)
+ E->viscosity.cdepv_ff[i] = 1.0; /* flavor factors for CDEPV */
+
/* read in information */
input_boolean("VISC_UPDATE",&(E->viscosity.update_allowed),"on",m);
input_int("rheol",&(E->viscosity.RHEOL),"3",m);
@@ -96,6 +99,21 @@
input_float("sdepv_misfit",&(E->viscosity.sdepv_misfit),"0.001",m);
+ 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->trace.nflavors > 10)
+ myerror(E,"error: too many flavors for CDEPV");
+ /* read in flavor factors */
+ input_float_vector("cdepv_ff",E->trace.nflavors,
+ (E->viscosity.cdepv_ff),m);
+ /* and take the log because we're using a geometric avg */
+ for(i=0;i<E->trace.nflavors;i++)
+ E->viscosity.cdepv_ff[i] = log(E->viscosity.cdepv_ff[i]);
+ }
+
+
input_boolean("low_visc_channel",&(E->viscosity.channel),"off",m);
input_boolean("low_visc_wedge",&(E->viscosity.wedge),"off",m);
@@ -146,13 +164,17 @@
void visc_from_mat();
void visc_from_T();
void visc_from_S();
+
+ void visc_from_P();
+ void visc_from_C();
+
void apply_viscosity_smoother();
void visc_to_node_interpolate();
void visc_from_nodes_to_gint();
void visc_from_gint_to_nodes();
- void visc_from_P();
+
int i,j,m;
float temp1,temp2,*vvvis;
double *TG;
@@ -164,13 +186,22 @@
else
visc_from_mat(E,evisc);
+ if(E->viscosity.CDEPV) /* compositional prefactor */
+ visc_from_C(E,evisc);
+
if(E->viscosity.SDEPV)
- visc_from_S(E,evisc,propogate);
-
+ visc_from_S(E,evisc,propogate);
+
if(E->viscosity.PDEPV) /* "plasticity" */
visc_from_P(E,evisc);
+ /* i think this should me placed differently i.e. before the
+ stress dependence but I won't change it because it's by
+ someone else
+
+ TWB
+ */
if(E->viscosity.channel || E->viscosity.wedge)
apply_low_visc_wedge_channel(E, evisc);
@@ -617,8 +648,53 @@
return;
}
+/*
+multiply with compositional factor which is determined by a geometric
+mean average from the tracer composition, assuming two flavors and
+compositions between zero and unity
+*/
+void visc_from_C( E, EEta)
+ struct All_variables *E;
+ float **EEta;
+{
+ float comp,comp_fac,CC[9],tcomp;
+ double vmean,cc_loc;
+ int m,l,z,jj,kk,i;
+
+ const int vpts = vpoints[E->mesh.nsd];
+ const int nel = E->lmesh.nel;
+ const int ends = enodes[E->mesh.nsd];
+ if(E->trace.nflavors != 2)
+ myerror(E,"sorry, CDEPV only supports two flavors");
+
+ for(m=1;m <= E->sphere.caps_per_proc;m++) {
+
+ for(i = 1; i <= nel; i++){
+ /* determine composition of each of the nodes of the
+ element */
+ for(kk = 1; kk <= ends; kk++){
+ CC[kk] = E->composition.comp_node[m][E->ien[m][i].node[kk]];
+ if(CC[kk] < 0)CC[kk]=0.0;
+ if(CC[kk] > 1)CC[kk]=1.0;
+ }
+ for(jj = 1; jj <= vpts; jj++){
+ /* compute mean composition */
+ cc_loc = 0.0;
+ for(kk = 1; kk <= ends; kk++)
+ cc_loc += CC[kk] * E->N.vpt[GNVINDEX(kk, jj)];
+
+ /* geometric mean of viscosity */
+ vmean = exp(cc_loc * E->viscosity.cdepv_ff[1] +
+ (1.0-cc_loc) * E->viscosity.cdepv_ff[0]);
+ /* multiply the viscosity with this prefactor */
+ EEta[m][ (i-1)*vpts + jj ] *= vmean;
+ } /* end jj loop */
+ } /* end el loop */
+ } /* end cap */
+}
+
void strain_rate_2_inv(E,m,EEDOT,SQRT)
struct All_variables *E;
float *EEDOT;
Modified: mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2007-08-09 00:36:01 UTC (rev 7789)
+++ mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2007-08-09 01:27:32 UTC (rev 7790)
@@ -86,6 +86,8 @@
float sdepv_trns[40];
+ int CDEPV; /* compositional viscosity */
+ float cdepv_ff[10]; /* flavor factors */
int PDEPV; /* "plasticity" law parameters */
float pdepv_a[40], pdepv_b[40], pdepv_y[40],pdepv_offset;
More information about the cig-commits
mailing list