[cig-commits] r7785 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Tue Aug 7 22:58:49 PDT 2007
Author: becker
Date: 2007-08-07 22:58:48 -0700 (Tue, 07 Aug 2007)
New Revision: 7785
Modified:
mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
mc/3D/CitcomS/trunk/lib/Drive_solvers.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Makefile.am
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/lib/output.h
mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
Log:
- added functions rtp2xyz, calc_cbase_at_tp, and
convert_pvec_to_cvec for spherical to Cartesian vector
conversion to Pan_problem_functions.c
Added "safe_malloc" to the same file because I didn't like the way
Malloc1 uses int as arguments (64 bit issues).
- added support for gzipped ascii output into subdirectories
of data_dir, option gzdir. To use, need to compile with
USE_GZDIR and link with the -lz zlib libraries.
Within gzdir output mode, selected by
output_format=gzdir
user can define vtk I/O (gzdir_vtkio=1) to write Cartesian
coordinates and velocities that are easy to combine to a single vtk
file in a post-processing step.
If gzlib_vtkio=0 the gzdir output option is meant to be identical to
the original ascii output, but files are all gzipped on the fly, and
the output that changes with each cycle is placed into
subdirectories of data_dir/ called data_dir/cycles/
For gzdir I/O, data_name is not used at all, all single output files
are placed in the data_dir/ directory without prefix.
- added a new rheology option, PDEPV, for pseudo-plasticity and
added a corresponding description to the manual.
- Modified Drive_solver stderr output to show v.v, and not only dv.dv/v.v
Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-08-08 03:37:43 UTC (rev 7784)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-08-08 05:58:48 UTC (rev 7785)
@@ -151,7 +151,7 @@
if(E->control.filter_temperature)
filter(E);
-
+
/* time0= CPU_time0()-time1; */
/* if(E->control.verbose) */
/* fprintf(E->fp_out,"time=%f\n",time0); */
Modified: mc/3D/CitcomS/trunk/lib/Drive_solvers.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-08-08 03:37:43 UTC (rev 7784)
+++ mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-08-08 05:58:48 UTC (rev 7785)
@@ -95,7 +95,7 @@
solve_constrained_flow_iterative(E);
- if (E->viscosity.SDEPV) {
+ if (E->viscosity.SDEPV || E->viscosity.PDEPV) {
for (m=1;m<=E->sphere.caps_per_proc;m++) {
delta_U[m] = (float *)malloc((neq+2)*sizeof(float));
@@ -118,13 +118,15 @@
Udot_mag = sqrt(global_fvdot(E,oldU,oldU,E->mesh.levmax));
dUdot_mag = vnorm_nonnewt(E,delta_U,oldU,E->mesh.levmax);
+
if(E->parallel.me==0){
- fprintf(stderr,"Stress dependent viscosity: DUdot = %.4e (%.4e) for iteration %d\n",dUdot_mag,Udot_mag,count);
- fprintf(E->fp,"Stress dependent viscosity: DUdot = %.4e (%.4e) for iteration %d\n",dUdot_mag,Udot_mag,count);
+ fprintf(stderr,"Stress dep. visc./plast.: DUdot = %.4e (%.4e) for iteration %d\n",
+ dUdot_mag,Udot_mag,count);
+ fprintf(E->fp,"Stress dep. visc./plast.: DUdot = %.4e (%.4e) for iteration %d\n",
+ 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]);
@@ -140,7 +142,7 @@
free((void *) delta_U[m]);
}
- } /*end if SDEPV*/
+ } /*end if SDEPV or PDEPV */
return;
}
@@ -188,7 +190,7 @@
}
solve_constrained_flow_iterative_pseudo_surf(E);
- if (E->viscosity.SDEPV) {
+ if (E->viscosity.SDEPV || E->viscosity.PDEPV) {
for (m=1;m<=E->sphere.caps_per_proc;m++) {
delta_U[m] = (float *)malloc((neq+2)*sizeof(float));
@@ -217,7 +219,7 @@
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]);
@@ -232,7 +234,7 @@
free((void *) delta_U[m]);
}
- } /*end if SDEPV */
+ } /*end if SDEPV or PDEPV */
E->monitor.topo_loop++;
}
get_STD_freesurf(E,E->slice.freesurf);
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-08 03:37:43 UTC (rev 7784)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-08 05:58:48 UTC (rev 7785)
@@ -944,7 +944,10 @@
char logfile[255];
E->fp = NULL;
- sprintf(logfile,"%s.log", E->control.data_file);
+ if (strcmp(E->output.format, "gzdir") == 0)
+ sprintf(logfile,"%s/log", E->control.data_dir);
+ else
+ sprintf(logfile,"%s.log", E->control.data_file);
E->fp = output_open(logfile);
return;
@@ -957,6 +960,9 @@
E->fptime = NULL;
if (E->parallel.me == 0) {
+ if (strcmp(E->output.format, "gzdir") == 0)
+ sprintf(timeoutput,"%s/time", E->control.data_dir);
+ else
sprintf(timeoutput,"%s.time", E->control.data_file);
E->fptime = output_open(timeoutput);
}
@@ -971,7 +977,10 @@
E->fp_out = NULL;
if (E->control.verbose) {
- sprintf(output_file,"%s.info.%d", E->control.data_file, E->parallel.me);
+ if (strcmp(E->output.format, "gzdir") == 0)
+ sprintf(output_file,"%s/info.%d", E->control.data_dir, E->parallel.me);
+ else
+ sprintf(output_file,"%s.info.%d", E->control.data_file, E->parallel.me);
E->fp_out = output_open(output_file);
}
@@ -1182,14 +1191,27 @@
}
else if (strcmp(E->output.format, "hdf5") == 0)
E->problem_output = h5output;
+#ifdef USE_GZDIR
+ else if (strcmp(E->output.format, "gzdir") == 0)
+ E->problem_output = gzdir_output;
else {
/* indicate error here */
if (E->parallel.me == 0) {
- fprintf(stderr, "wrong output_format, must be either 'ascii' or 'hdf5'\n");
- fprintf(E->fp, "wrong output_format, must be either 'ascii' or 'hdf5'\n");
+ fprintf(stderr, "wrong output_format, must be 'ascii', 'hdf5', or 'gzdir'\n");
+ fprintf(E->fp, "wrong output_format, must be 'ascii', 'hdf5', or 'gzdir'\n");
}
parallel_process_termination(E);
}
+#else
+ else {
+ /* indicate error here */
+ if (E->parallel.me == 0) {
+ fprintf(stderr, "wrong output_format, must be 'ascii' or 'gzdir' (USE_GZDIR undefined)\n");
+ fprintf(E->fp, "wrong output_format, must be 'ascii' or 'gzdir' (USE_GZDIR undefined)\n");
+ }
+ parallel_process_termination(E);
+ }
+#endif
output_parse_optional(E);
}
Modified: mc/3D/CitcomS/trunk/lib/Makefile.am
===================================================================
--- mc/3D/CitcomS/trunk/lib/Makefile.am 2007-08-08 03:37:43 UTC (rev 7784)
+++ mc/3D/CitcomS/trunk/lib/Makefile.am 2007-08-08 05:58:48 UTC (rev 7785)
@@ -85,6 +85,7 @@
Nodal_mesh.c \
Output.c \
output.h \
+ Output_gzdir.c \
Output_h5.c \
output_h5.h \
Pan_problem_misc_functions.c \
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2007-08-08 03:37:43 UTC (rev 7784)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2007-08-08 05:58:48 UTC (rev 7785)
@@ -64,6 +64,14 @@
input_string("output_format", E->output.format, "ascii",m);
input_string("output_optional", E->output.optional, "surf,botm,tracer,comp_el",m);
+ /* gzdir type of I/O */
+ if(strcmp(E->output.format, "gzdir") == 0){
+ input_boolean("gzdir_vtkio",&(E->output.gzdir_vtkio),"off",m);
+ E->output.gzdir_vtkbase_init = 0;
+ E->output.gzdir_vtkbase_save = 1; /* should we save the basis vectors? (memory!) */
+ fprintf(stderr,"gzdir: vtkio: %i save basis vectors: %i\n",
+ E->output.gzdir_vtkio,E->output.gzdir_vtkbase_save);
+ }
}
@@ -76,6 +84,7 @@
/*output_mat(E);*/
}
+
output_velo(E, cycles);
output_visc(E, cycles);
Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-08 03:37:43 UTC (rev 7784)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-08 05:58:48 UTC (rev 7785)
@@ -51,7 +51,13 @@
#include "phase_change.h"
#include "parallel_related.h"
+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 get_process_identifier()
{
int pid;
@@ -385,3 +391,72 @@
{
return 1.0;
}
+
+/* convert r,theta,phi system to cartesian, xout[3] */
+void rtp2xyz(float r, float theta, float phi, float *xout)
+{
+ float rst;
+ rst = r * sin(theta);
+ xout[0] = rst * cos(phi); /* x */
+ xout[1] = rst * sin(phi); /* y */
+ xout[2] = r * cos(theta);
+}
+
+
+/* compute base vectors for conversion of polar to cartesian vectors
+ base[9]
+*/
+void calc_cbase_at_tp(float theta, float phi, float *base)
+{
+
+
+ double ct,cp,st,sp;
+
+ ct=cos(theta);
+ cp=cos(phi);
+ st=sin(theta);
+ sp=sin(phi);
+ /* r */
+ base[0]= st * cp;
+ base[1]= st * sp;
+ base[2]= ct;
+ /* theta */
+ base[3]= ct * cp;
+ base[4]= ct * sp;
+ base[5]= -st;
+ /* phi */
+ base[6]= -sp;
+ base[7]= cp;
+ base[8]= 0.0;
+}
+
+/* given a base from calc_cbase_at_tp, convert a polar vector to
+ cartesian */
+void convert_pvec_to_cvec(float vr,float vt,
+ float vp, float *base,
+ float *cvec)
+{
+ int i;
+ for(i=0;i<3;i++){
+ cvec[i] = base[i] * vr;
+ cvec[i] += base[3+i]* vt;
+ cvec[i] += base[6+i]* vp;
+ }
+}
+/*
+ like malloc, but with test
+
+ similar to Malloc1 but I didn't like the int as argument
+
+*/
+void *safe_malloc (size_t size)
+{
+ void *tmp;
+
+ if ((tmp = malloc(size)) == NULL) {
+ fprintf(stderr, "safe_malloc: could not allocate memory, %.3f MB\n",
+ (float)size/(1024*1024.));
+ parallel_process_termination();
+ }
+ return (tmp);
+}
Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-08-08 03:37:43 UTC (rev 7784)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-08-08 05:58:48 UTC (rev 7785)
@@ -239,9 +239,9 @@
count, CPU_time0()-time0, 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 dv/v=%.3e"
+ fprintf(stderr,"AhatP (%03d) after %.3f seconds with 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, E->monitor.vdotv ,E->monitor.incompressibility,
0.0, 0.0, E->monitor.solution_cycles);
}
@@ -311,9 +311,9 @@
" and dp/p=%.3e for step %d\n",
count, CPU_time0()-time0, E->monitor.incompressibility,
dvelocity, dpressure, E->monitor.solution_cycles);
- fprintf(stderr, "AhatP (%03d) after %g seconds with div/v=%.3e dv/v=%.3e"
+ fprintf(stderr, "AhatP (%03d) after %.3f seconds with 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, E->monitor.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-08-08 03:37:43 UTC (rev 7784)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-08 05:58:48 UTC (rev 7785)
@@ -52,6 +52,12 @@
E->viscosity.T[i] = 0.0;
E->viscosity.Z[i] = 0.0;
E->viscosity.E[i] = 0.0;
+
+ E->viscosity.pdepv_a[i] = 1.e20; /* \sigma_y = min(a + b * (1-r),y) */
+ E->viscosity.pdepv_b[i] = 0.0;
+ E->viscosity.pdepv_y[i] = 1.e20;
+
+
}
/* read in information */
@@ -71,10 +77,25 @@
E->viscosity.sdepv_misfit = 1.0;
input_boolean("SDEPV",&(E->viscosity.SDEPV),"off",m);
if (E->viscosity.SDEPV) {
- input_float("sdepv_misfit",&(E->viscosity.sdepv_misfit),"0.001",m);
- input_float_vector("sdepv_expt",E->viscosity.num_mat,(E->viscosity.sdepv_expt),m);
+ input_float_vector("sdepv_expt",E->viscosity.num_mat,(E->viscosity.sdepv_expt),m);
}
+
+ input_boolean("PDEPV",&(E->viscosity.PDEPV),"off",m); /* plasticity addition by TWB */
+ if (E->viscosity.PDEPV) {
+ E->viscosity.pdepv_visited = 0;
+
+ input_boolean("pdepv_eff",&(E->viscosity.pdepv_eff),"on",m);
+ input_float_vector("pdepv_a",E->viscosity.num_mat,(E->viscosity.pdepv_a),m);
+ input_float_vector("pdepv_b",E->viscosity.num_mat,(E->viscosity.pdepv_b),m);
+ input_float_vector("pdepv_y",E->viscosity.num_mat,(E->viscosity.pdepv_y),m);
+
+ input_float("pdepv_offset",&(E->viscosity.pdepv_offset),"0.0",m);
+ }
+ if(E->viscosity.PDEPV || E->viscosity.SDEPV)
+ input_float("sdepv_misfit",&(E->viscosity.sdepv_misfit),"0.001",m);
+
+
input_boolean("low_visc_channel",&(E->viscosity.channel),"off",m);
input_boolean("low_visc_wedge",&(E->viscosity.wedge),"off",m);
@@ -130,6 +151,7 @@
void visc_from_nodes_to_gint();
void visc_from_gint_to_nodes();
+ void visc_from_P();
int i,j,m;
float temp1,temp2,*vvvis;
@@ -145,7 +167,13 @@
if(E->viscosity.SDEPV)
visc_from_S(E,evisc,propogate);
+ if(E->viscosity.PDEPV){ /* "plasticity" */
+ //visc_from_P(E,evisc);
+ if(E->parallel.me == 0)
+ fprintf(stderr,"PDEPV switched off to test the stress dependent loop\n");
+ }
+
if(E->viscosity.channel || E->viscosity.wedge)
apply_low_visc_wedge_channel(E, evisc);
@@ -457,7 +485,7 @@
two = 2.0;
for(m=1;m<=E->sphere.caps_per_proc;m++) {
- strain_rate_2_inv(E,m,eedot,1);
+ strain_rate_2_inv(E,m,eedot,1); /* should there be a check if velocities have been computed here? */
for(e=1;e<=nel;e++) {
exponent1= one/E->viscosity.sdepv_expt[E->mat[m][e]-1];
@@ -471,8 +499,105 @@
return;
}
+void visc_from_P(E,EEta) /* "plasticity" implementation
+
+ viscosity will be limited by a yield stress
+
+ \sigma_y = min(a + b * (1-r), y)
+
+ where a,b,y are parameters input via pdepv_a,b,y
+
+ and
+
+ \eta_y = \sigma_y / (2 \eps_II)
+
+ where \eps_II is the second invariant. Then
+
+ \eta_eff = (\eta_0 \eta_y)/(\eta_0 + \eta_y)
+
+ for pdepv_eff = 1
+
+ or
+ \eta_eff = min(\eta_0,\eta_y)
+
+ for pdepv_eff = 0
+
+ where \eta_0 is the regular viscosity
+
+
+ TWB
+
+ */
+ struct All_variables *E;
+ float **EEta;
+{
+ float *eedot,zz[9],zzz,tau,eta_p,eta_new;
+ int m,e,l,z,jj,kk;
+ const int vpts = vpoints[E->mesh.nsd];
+ const int nel = E->lmesh.nel;
+ const int ends = enodes[E->mesh.nsd];
+
+ void strain_rate_2_inv();
+
+ eedot = (float *) malloc((2+nel)*sizeof(float));
+
+ for(m=1;m<=E->sphere.caps_per_proc;m++) {
+
+ if(E->viscosity.pdepv_visited){
+
+ strain_rate_2_inv(E,m,eedot,1); /* get second invariant for all elements */
+
+ }else{
+ for(e=1;e<=nel;e++) /* initialize with unity if no velocities around */
+ eedot[e] = 1.0e-5;
+ if(m == E->sphere.caps_per_proc)
+ E->viscosity.pdepv_visited = 1;
+ }
+
+ for(e=1;e<=nel;e++) { /* loop through all elements */
+
+ l = E->mat[m][e]; /* material of this element */
+
+ for(kk=1;kk <= ends;kk++) /* nodal depths */
+ zz[kk] = (1.0 - E->sx[m][3][E->ien[m][e].node[kk]]); /* for depth, zz = 1 - r */
+
+ for(jj=1;jj <= vpts;jj++){ /* loop through integration points */
+
+ zzz = 0.0; /* get mean depth of integration point */
+ for(kk=1;kk<=ends;kk++)
+ zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+
+ /* depth dependent yield stress */
+ tau = E->viscosity.pdepv_a[l] + zzz * E->viscosity.pdepv_b[l];
+
+ /* min of depth dep. and constant yield stress */
+ tau = min(tau, E->viscosity.pdepv_y[l]);
+
+ /* yield viscosity */
+ eta_p = tau/(2.0 * eedot[e] + 1e-7) + E->viscosity.pdepv_offset;
+
+
+ if(E->viscosity.pdepv_eff){
+ /* two dashpots in series */
+ eta_new = 1.0/(1.0/EEta[m][ (e-1)*vpts + jj ] + 1.0/eta_p);
+ }else{
+ /* min viscosities*/
+ eta_new = min(EEta[m][ (e-1)*vpts + jj ], eta_p);
+ }
+ //fprintf(stderr,"%11g %11g %11g %11g %11g\n",eedot[e],tau,eta_p,eta_new,EEta[m][(e-1)*vpts + jj]);
+ EEta[m][(e-1)*vpts + jj] = eta_new;
+ } /* end integration point loop */
+ } /* end element loop */
+
+ } /* end caps loop */
+ free ((void *)eedot);
+ return;
+}
+
+
+
void strain_rate_2_inv(E,m,EEDOT,SQRT)
struct All_variables *E;
float *EEDOT;
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-08 03:37:43 UTC (rev 7784)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-08 05:58:48 UTC (rev 7785)
@@ -658,6 +658,12 @@
int comp_el; /* whether to output composition at elements */
int comp_nd; /* whether to output composition at nodes */
+
+#ifdef USE_GZDIR
+ /* flags only used by GZDIR */
+ int gzdir_vtkio,gzdir_vtkbase_init,gzdir_vtkbase_save;
+ float *gzdir_vtkbase;
+#endif
};
Modified: mc/3D/CitcomS/trunk/lib/output.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/output.h 2007-08-08 03:37:43 UTC (rev 7784)
+++ mc/3D/CitcomS/trunk/lib/output.h 2007-08-08 05:58:48 UTC (rev 7785)
@@ -43,4 +43,7 @@
}
#endif
+#ifdef USE_GZDIR
+void gzdir_output(struct All_variables *, int );
#endif
+#endif
Modified: mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2007-08-08 03:37:43 UTC (rev 7784)
+++ mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2007-08-08 05:58:48 UTC (rev 7785)
@@ -85,6 +85,12 @@
float sdepv_expt[40];
float sdepv_trns[40];
+
+
+ int PDEPV; /* "plasticity" law parameters */
+ float pdepv_a[40], pdepv_b[40], pdepv_y[40],pdepv_offset;
+ int pdepv_eff,pdepv_visited;
+
int TDEPV;
int TDEPV_AVE;
float N0[40];
More information about the cig-commits
mailing list