[cig-commits] r7788 - in mc/3D/CitcomS/trunk: doc/manual lib
becker at geodynamics.org
becker at geodynamics.org
Wed Aug 8 15:55:40 PDT 2007
Author: becker
Date: 2007-08-08 15:55:39 -0700 (Wed, 08 Aug 2007)
New Revision: 7788
Added:
mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
Modified:
mc/3D/CitcomS/trunk/doc/manual/citcoms.lyx
mc/3D/CitcomS/trunk/lib/Composition_related.c
mc/3D/CitcomS/trunk/lib/Construct_arrays.c
mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
mc/3D/CitcomS/trunk/lib/General_matrix_functions.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/Output_gzdir.c
mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/trunk/lib/Tracer_setup.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/lib/tracer_defs.h
Log:
- added new routine Ggrd_handling which will take care of netcdf grd
input. For now, only the tracer innitial conditions can be set by
using, for example,
ictracer_grd_layers=2
ictracer_grd_file="tmask.grd"
for this to work, the ggrd libraries from HC, as well as the -lgmt
and -lnetcdf need to be linked, and -DUSE_GGRD needs to be defined
on compile time.
- rewrote TDEPV viscosity options 3 and 4 slightly to make them more
efficient
- added TDEPV viscosity option 6
eta = N_0 exp(E(T_0-T) + (1-z) Z_0 )
- took layers subroutine in Construct_arrays.c apart to allow calling
the same function in a more flexible from other subroutines
- addded error handling routine "myerror" to Pan_problem_misc_functions
- made the output behavior of the multigrid solver convergence
parameters consisent with verbosity setting (ie. use report instead
of always stderr output)
- changed tracer input parameter "tracer" and "chemical_buoyancy" to
boolean routine for added flexibility (I hope)
- made check_initial_composition consistent with tracer advection,
exits after warning only for E->trace.itracer_warnings==1
- debugging output in Full_tracer_advection.c and Tracer_setup only
gets printed to file for E->control.verbose == 1
- fixed some typos in Output_gzdir (that might have been introduced by
asynchronous svn updating or something like that)
Modified: mc/3D/CitcomS/trunk/doc/manual/citcoms.lyx
===================================================================
--- mc/3D/CitcomS/trunk/doc/manual/citcoms.lyx 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/doc/manual/citcoms.lyx 2007-08-08 22:55:39 UTC (rev 7788)
@@ -1,4 +1,4 @@
-#LyX 1.4.3 created this file. For more info see http://www.lyx.org/
+#LyX 1.4.4 created this file. For more info see http://www.lyx.org/
\lyxformat 245
\begin_document
\begin_header
@@ -11077,7 +11077,22 @@
\end_layout
+\begin_layout Standard
+When
+\family typewriter
+rheol
+\family default
+=6, temperature-dependent viscosity is computed by
+\family typewriter
+\size footnotesize
+
+\begin_inset Formula $\eta=\eta_{0}\times\exp(\eta_{E}(\eta_{T}-T)+(1-z)\eta_{Z})$
\end_inset
+
+
+\end_layout
+
+\end_inset
</cell>
</row>
<row topline="true">
Modified: mc/3D/CitcomS/trunk/lib/Composition_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Composition_related.c 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/Composition_related.c 2007-08-08 22:55:39 UTC (rev 7788)
@@ -45,10 +45,11 @@
{
int m = E->parallel.me;
- input_int("chemical_buoyancy",&(E->composition.ichemical_buoyancy),
- "1,0,nomax",m);
+ input_boolean("chemical_buoyancy",
+ &(E->composition.ichemical_buoyancy),
+ "1,0,nomax",m);
- if (E->composition.ichemical_buoyancy==1) {
+ if (E->composition.ichemical_buoyancy) {
input_double("buoyancy_ratio",
&(E->composition.buoyancy_ratio),"1.0",m);
@@ -106,7 +107,7 @@
{
E->composition.on = 0;
- if (E->composition.ichemical_buoyancy==1 ||
+ if (E->composition.ichemical_buoyancy ||
E->composition.icompositional_rheology)
E->composition.on = 1;
@@ -117,16 +118,17 @@
parallel_process_termination();
}
- if (E->composition.ichemical_buoyancy==0)
- fprintf(E->trace.fpt,"Passive Tracers\n");
+ if (!E->composition.ichemical_buoyancy)
+ fprintf(E->trace.fpt,"Passive Tracers\n");
+ else
+ fprintf(E->trace.fpt,"Active Tracers\n");
- if (E->composition.ichemical_buoyancy==1)
- fprintf(E->trace.fpt,"Active Tracers\n");
+ if (E->composition.ibuoy_type==1)
+ fprintf(E->trace.fpt,"Ratio Method\n");
+ if (E->composition.ibuoy_type==0)
+ fprintf(E->trace.fpt,"Absolute Method\n");
- if (E->composition.ibuoy_type==1) fprintf(E->trace.fpt,"Ratio Method\n");
- if (E->composition.ibuoy_type==0) fprintf(E->trace.fpt,"Absolute Method\n");
-
fprintf(E->trace.fpt,"Buoyancy Ratio: %f\n", E->composition.buoyancy_ratio);
if (E->composition.ireset_initial_composition==0)
@@ -209,12 +211,13 @@
static void init_composition(struct All_variables *E)
{
- if (E->composition.ichemical_buoyancy==1 && E->composition.ibuoy_type==1) {
- fill_composition(E);
- check_initial_composition(E);
- init_bulk_composition(E);
- }
- return;
+ if (E->composition.ichemical_buoyancy &&
+ E->composition.ibuoy_type) {
+ fill_composition(E);
+ check_initial_composition(E);
+ init_bulk_composition(E);
+ }
+ return;
}
@@ -226,7 +229,7 @@
fprintf(E->trace.fpt,"WARNING(check_initial_composition)-number of tracers is REALLY LOW\n");
fflush(E->trace.fpt);
fprintf(stderr,"WARNING(check_initial_composition)-number of tracers is REALLY LOW\n");
- exit(10);
+ if (E->trace.itracer_warnings==1) exit(10); /* made this consistent with tracer advection */
}
}
Modified: mc/3D/CitcomS/trunk/lib/Construct_arrays.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Construct_arrays.c 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/Construct_arrays.c 2007-08-08 22:55:39 UTC (rev 7788)
@@ -30,6 +30,8 @@
#include "element_definitions.h"
#include "global_defs.h"
+int layers_r(struct All_variables *,float );
+
/*========================================================
Function to make the IEN array for a mesh of given
dimension. IEN is an externally defined structure array
@@ -746,24 +748,32 @@
int layers(E,m,node)
- struct All_variables *E;
- int m,node;
+ struct All_variables *E;
+ int m,node;
{
- float zlith, z410, zlm;
-
- int llayers = 0;
- zlith=E->viscosity.zlith;
- z410=E->viscosity.z410;
- zlm=E->viscosity.zlm;
-
- if (E->sx[m][3][node]>(E->sphere.ro-zlith))
- llayers = 1;
- else if ((E->sx[m][3][node])>(E->sphere.ro-z410) && E->sx[m][3][node]<=(E->sphere.ro-zlith))
- llayers = 2;
- else if ((E->sx[m][3][node])>(E->sphere.ro-zlm) && E->sx[m][3][node]<=(E->sphere.ro-z410))
- llayers = 3;
- else
- llayers = 4;
-
- return (llayers);
- }
+ return(layers_r(E,E->sx[m][3][node]));
+}
+/* took this apart to allow call from other subroutines */
+int layers_r(E,r)
+ struct All_variables *E;
+ float r;
+{
+ float zlith, z410, zlm;
+
+ int llayers = 0;
+ zlith = E->viscosity.zlith;
+ z410 = E->viscosity.z410;
+ zlm = E->viscosity.zlm;
+
+ if (r > (E->sphere.ro-zlith))
+ llayers = 1;
+ else if ((r > (E->sphere.ro-z410))&&
+ (r <= (E->sphere.ro-zlith)))
+ llayers = 2;
+ else if ((r > (E->sphere.ro-zlm)) &&
+ (r <= (E->sphere.ro-z410)))
+ llayers = 3;
+ else
+ llayers = 4;
+ return (llayers);
+}
Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-08-08 22:55:39 UTC (rev 7788)
@@ -305,23 +305,25 @@
parallel_process_sync(E);
- fprintf(E->trace.fpt, "Entering lost_souls()\n");
+ if(E->control.verbose)
+ fprintf(E->trace.fpt, "Entering lost_souls()\n");
E->trace.istat_isend=E->trace.ilater[j];
-
-
/* debug */
- for (kk=1; kk<=E->trace.istat_isend; kk++) {
+ if(E->control.verbose){
+ for (kk=1; kk<=E->trace.istat_isend; kk++) {
fprintf(E->trace.fpt, "tracer#=%d xx=(%g,%g,%g)\n", kk,
E->trace.rlater[j][0][kk],
E->trace.rlater[j][1][kk],
E->trace.rlater[j][2][kk]);
+ }
+ fflush(E->trace.fpt);
}
- fflush(E->trace.fpt);
/**/
+
/* initialize isend and ireceive */
/* # of neighbors in the horizontal plane */
isize[j]=E->trace.ilater[j]*E->trace.number_of_tracer_quantities;
@@ -388,16 +390,18 @@
/** debug **/
- for (kk=0;kk<=num_ngb;kk++) {
+ if(E->control.verbose){
+ for (kk=0;kk<=num_ngb;kk++) {
if(kk==0)
- isource_proc=E->parallel.me;
+ isource_proc=E->parallel.me;
else
- isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
-
- fprintf(E->trace.fpt,"%d send %d to proc %d\n",
- E->parallel.me,isend[j][kk],isource_proc);
- fprintf(E->trace.fpt,"%d recv %d from proc %d\n",
- E->parallel.me,ireceive[j][kk],isource_proc);
+ isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+
+ fprintf(E->trace.fpt,"%d send %d to proc %d\n",
+ E->parallel.me,isend[j][kk],isource_proc);
+ fprintf(E->trace.fpt,"%d recv %d from proc %d\n",
+ E->parallel.me,ireceive[j][kk],isource_proc);
+ }
}
/**/
@@ -745,9 +749,10 @@
E->trace.ielement[j][E->trace.ntracers[j]]=iel;
}
-
- fprintf(E->trace.fpt,"Freeing memory in lost_souls()\n");
- fflush(E->trace.fpt);
+ if(E->control.verbose){
+ fprintf(E->trace.fpt,"Freeing memory in lost_souls()\n");
+ fflush(E->trace.fpt);
+ }
parallel_process_sync(E);
/* Free Arrays */
@@ -759,8 +764,10 @@
free(receive[j][kk]);
}
- fprintf(E->trace.fpt,"Leaving lost_souls()\n");
- fflush(E->trace.fpt);
+ if(E->control.verbose){
+ fprintf(E->trace.fpt,"Leaving lost_souls()\n");
+ fflush(E->trace.fpt);
+ }
return;
}
@@ -1865,6 +1872,13 @@
fprintf(E->trace.fpt,"Layered tracer flavors\n");
fprintf(E->trace.fpt,"Interface Height: %f\n",E->trace.z_interface);
}
+#ifdef USE_GGRD
+ else if(E->trace.ic_method_for_flavors == 1) {
+ fprintf(E->trace.fpt,"netcdf grd assigned tracer flavors\n");
+ fprintf(E->trace.fpt,"file: %s top %i layeres\n",E->trace.ggrd_file,
+ E->trace.ggrd_layers);
+ }
+#endif
else {
fprintf(E->trace.fpt,"Sorry-This IC methods for Flavors are Unavailable %d\n",E->trace.ic_method_for_flavors);
fflush(E->trace.fpt);
Modified: mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/General_matrix_functions.c 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/General_matrix_functions.c 2007-08-08 22:55:39 UTC (rev 7788)
@@ -208,6 +208,8 @@
int count,counts,cycles,convergent,valid;
int i, neq, gneq,m;
+ char message[200];
+
double CPU_time0(),initial_time,time;
double residual,prior_residual,r0;
double *D1[NCS], *r[NCS], *Au[NCS];
@@ -236,15 +238,21 @@
else {
counts =0;
- if(E->parallel.me==0) fprintf(stderr,"resi = %.6e for iter %d acc %.6e\n",residual,counts,acc);
- if(E->parallel.me==0) fprintf(E->fp,"resi = %.6e for iter %d acc %.6e\n",residual,counts,acc);
+ if(E->parallel.me==0){ /* output */
+ snprintf(message,200,"resi = %.6e for iter %d acc %.6e",residual,counts,acc);
+ record(E,message);
+ report(E,message);
+ }
do {
residual=multi_grid(E,d0,F,acc,high_lev);
valid = (residual < acc)?1:0;
counts ++;
- if(E->parallel.me==0) fprintf(stderr,"resi = %.6e for iter %d acc %.6e\n",residual,counts,acc);
- if(E->parallel.me==0) fprintf(E->fp,"resi = %.6e for iter %d acc %.6e\n",residual,counts,acc);
+ if(E->parallel.me==0){ /* output */
+ snprintf(message,200,"resi = %.6e for iter %d acc %.6e",residual,counts,acc);
+ record(E,message);
+ report(E,message);
+ }
} while (!valid);
cycles = counts;
Added: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2007-08-08 22:55:39 UTC (rev 7788)
@@ -0,0 +1,155 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *<LicenseText>
+ *
+ * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
+ * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
+ * Copyright (C) 1994-2005, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+/*
+
+routines that deal with GMT/netcdf grd I/O as supported through
+the ggrd subroutines of the hc package
+
+*/
+
+#include <math.h>
+#include "global_defs.h"
+#include "parsing.h"
+#include "parallel_related.h"
+#include "composition_related.h"
+
+#ifdef USE_GGRD
+
+#include "hc.h" /* ggrd and hc packages */
+
+void ggrd_init_tracer_flavors(struct All_variables *);
+int layers_r(struct All_variables *,float );
+void myerror(struct All_variables *,char *);
+
+/*
+
+assign tracer flavor based on its depth (within top n layers),
+and the grd value
+
+
+*/
+void ggrd_init_tracer_flavors(struct All_variables *E)
+{
+ int j, kk, number_of_tracers;
+ double rad,theta,phi,indbl;
+ char char_dummy[1],error[255],gmt_bc[10];
+ struct ggrd_gt ggrd_ict[1];
+ /* for dealing with several processors */
+ MPI_Status mpi_stat;
+ int mpi_rc;
+ int mpi_inmsg, mpi_success_message = 1;
+
+ report(E,"ggrd_init_tracer_flavors: ggrd mat init");
+
+ /*
+ are we global?
+ */
+ if (E->parallel.nprocxy == 12){
+ /* use GMT's geographic boundary conditions */
+ sprintf(gmt_bc,"-Lg");
+ }else{ /* regional */
+ sprintf(gmt_bc,"");
+ }
+
+ /*
+ initialize the ggrd control
+ */
+ if(E->parallel.me > 0){
+ /* wait for previous processor */
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1),
+ 0, MPI_COMM_WORLD, &mpi_stat);
+ }
+ if(ggrd_grdtrack_init_general(FALSE,E->trace.ggrd_file,
+ char_dummy,gmt_bc,
+ ggrd_ict,FALSE,FALSE)){
+ myerror(E,"ggrd tracer init error");
+ }
+ if(E->parallel.me < E->parallel.nproc-1){
+ /* tell the next proc to go ahead */
+ mpi_rc = MPI_Send(&mpi_success_message, 1,
+ MPI_INT, (E->parallel.me+1), 0, MPI_COMM_WORLD);
+ }else{
+ report(E,"ggrd_init_tracer_flavors: last processor done with ggrd mat init");
+ }
+ /* init done */
+
+ /* assign values to each tracer based on grd file */
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ number_of_tracers = E->trace.ntracers[j];
+ for (kk=1;kk <= number_of_tracers;kk++) {
+ rad = E->trace.basicq[j][2][kk]; /* tracer radius */
+ if(layers_r(E,rad) <= E->trace.ggrd_layers){
+ /*
+ in top layers
+ */
+ phi = E->trace.basicq[j][1][kk];
+ theta = E->trace.basicq[j][0][kk];
+ /* interpolate from grid */
+ if(!ggrd_grdtrack_interpolate_tp((double)theta,(double)phi,
+ ggrd_ict,&indbl,FALSE)){
+ snprintf(error,255,"read_mat_from_ggrd_file: interpolation error at theta: %g phi: %g",
+ theta,phi);
+ myerror(E,error);
+ }
+ /* limit to 0 or 1 */
+ if(indbl < .5)
+ indbl = 0.0;
+ else
+ indbl = 1.0;
+ E->trace.extraq[j][0][kk]= indbl;
+ }else{
+ /* below */
+ E->trace.extraq[j][0][kk]=0.0;
+ }
+ }
+ }
+
+ /* free grd structure */
+ ggrd_grdtrack_free_gstruc(ggrd_ict);
+ report(E,"ggrd tracer init done");
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#endif
+
+
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-08 22:55:39 UTC (rev 7788)
@@ -631,7 +631,7 @@
for(i=0;i<=E->lmesh.neq;i++)
E->U[j][i] = E->temp[j][i] = E->temp1[j][i] = 0.0;
- if(E->control.tracer==1) {
+ if(E->control.tracer == 1) {
for(i=1;i<=E->mesh.nsd;i++) {
E->GV[j][i]=(float*) malloc(((E->lmesh.nno+1)*E->parallel.nproc+1)*sizeof(float));
E->GV1[j][i]=(float*) malloc(((E->lmesh.nno+1)*E->parallel.nproc+1)*sizeof(float));
Modified: mc/3D/CitcomS/trunk/lib/Makefile.am
===================================================================
--- mc/3D/CitcomS/trunk/lib/Makefile.am 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/Makefile.am 2007-08-08 22:55:39 UTC (rev 7788)
@@ -73,6 +73,7 @@
element_definitions.h \
General_matrix_functions.c \
global_defs.h \
+ Ggrd_handling.c \
Global_operations.c \
hdf5_related.h \
Initial_temperature.c \
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2007-08-08 22:55:39 UTC (rev 7788)
@@ -69,8 +69,8 @@
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);
+ //fprintf(stderr,"gzdir: vtkio: %i save basis vectors: %i\n",
+ // E->output.gzdir_vtkio,E->output.gzdir_vtkbase_save);
}
}
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-08-08 22:55:39 UTC (rev 7788)
@@ -464,7 +464,7 @@
char output_file[255];
gzFile *fp1;
- snprintf(output_file,255,"%s/%d/stress.%d.%d", E->control.data_dir,
+ snprintf(output_file,255,"%s/%d/stress.%d.%d.gz", E->control.data_dir,
cycles,E->parallel.me, cycles);
fp1 = gzdir_output_open(output_file,"w");
@@ -540,7 +540,7 @@
char output_file[255];
gzFile *fp1;
- snprintf(output_file,255,"%s/%d/pressure.%d.%d", E->control.data_dir,cycles,
+ snprintf(output_file,255,"%s/%d/pressure.%d.%d.gz", E->control.data_dir,cycles,
E->parallel.me, cycles);
fp1 = gzdir_output_open(output_file,"w");
@@ -565,7 +565,7 @@
char output_file[255];
gzFile *fp1;
- snprintf(output_file,255,"%s/%d/tracer.%d.%d", E->control.data_dir,
+ 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");
@@ -578,14 +578,14 @@
for(n=1;n<=E->trace.ntracers[j];n++) {
/* write basic quantities (coordinate) */
- gzprintf(fp1,"%.12e %.12e %.12e",
+ gzprintf(fp1,"%9.5e %9.5e %9.5e",
E->trace.basicq[j][0][n],
E->trace.basicq[j][1][n],
E->trace.basicq[j][2][n]);
/* write extra quantities */
for (i=0; i<E->trace.number_of_extra_quantities; i++) {
- gzprintf(fp1," %.12e", E->trace.extraq[j][i][n]);
+ gzprintf(fp1," %9.5e", E->trace.extraq[j][i][n]);
}
gzprintf(fp1, "\n");
}
@@ -603,7 +603,7 @@
char output_file[255];
gzFile *fp1;
- snprintf(output_file,255,"%s/%d/comp_nd.%d.%d",
+ snprintf(output_file,255,"%s/%d/comp_nd.%d.%d.gz",
E->control.data_dir,
cycles,
E->parallel.me, cycles);
@@ -633,7 +633,7 @@
char output_file[255];
gzFile *fp1;
- snprintf(output_file,255,"%s/%d/comp_el.%d.%d", E->control.data_dir,
+ snprintf(output_file,255,"%s/%d/comp_el.%d.%d.gz", E->control.data_dir,
cycles,E->parallel.me, cycles);
fp1 = gzdir_output_open(output_file,"w");
Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-08 22:55:39 UTC (rev 7788)
@@ -55,6 +55,7 @@
void rtp2xyz(float , float , float, float *);
void convert_pvec_to_cvec(float ,float , float , float *,float *);
void *safe_malloc (size_t );
+void myerror(struct All_variables *,char *);
@@ -147,7 +148,8 @@
buoy[m][i] = temp * E->T[m][i];
/* chemical buoyancy */
- if(E->control.tracer && E->composition.ichemical_buoyancy==1) {
+ if(E->control.tracer &&
+ (E->composition.ichemical_buoyancy)) {
temp2 = E->composition.buoyancy_ratio * temp;
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.nno;i++)
@@ -392,7 +394,11 @@
return 1.0;
}
-/* convert r,theta,phi system to cartesian, xout[3] */
+/* convert r,theta,phi system to cartesian, xout[3]
+ there's a double version of this in Tracer_setup called
+ sphere_to_cart
+
+*/
void rtp2xyz(float r, float theta, float phi, float *xout)
{
float rst;
@@ -460,3 +466,13 @@
}
return (tmp);
}
+/* error handling routine, TWB */
+
+void myerror(struct All_variables *E,char *message)
+{
+ E->control.verbose = 1;
+ record(E,message);
+ fprintf(stderr,"node %3i: error: %s\n",
+ E->parallel.me,message);
+ parallel_process_termination();
+}
Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-08-08 22:55:39 UTC (rev 7788)
@@ -45,6 +45,9 @@
#include "parallel_related.h"
#include "composition_related.h"
+#ifdef USE_GGRD
+void ggrd_init_tracer_flavors(struct All_variables *);
+#endif
int icheck_that_processor_shell(struct All_variables *E,
int j, int nprocessor, double rad);
@@ -77,7 +80,7 @@
void full_tracer_input();
int m=E->parallel.me;
- input_int("tracer",&(E->control.tracer),"0",m);
+ input_boolean("tracer",&(E->control.tracer),"0",m);
if(E->control.tracer) {
/* tracer_ic_method=0 (random generated array) */
@@ -85,8 +88,9 @@
/* tracer_ic_method=2 (each proc reads its restart file) */
input_int("tracer_ic_method",&(E->trace.ic_method),"0,0,nomax",m);
- if (E->trace.ic_method==0)
+ if (E->trace.ic_method==0){
input_int("tracers_per_element",&(E->trace.itperel),"10,0,nomax",m);
+ }
else if (E->trace.ic_method==1)
input_string("tracer_file",E->trace.tracer_file,"tracer.dat",m);
else if (E->trace.ic_method==2) {
@@ -108,9 +112,22 @@
input_int("tracer_flavors",&(E->trace.nflavors),"0,0,nomax",m);
- input_int("ic_method_for_flavors",&(E->trace.ic_method_for_flavors),"0,0,nomax",m);
- if (E->trace.ic_method_for_flavors == 0)
- input_double("z_interface",&(E->trace.z_interface),"0.7",m);
+ input_int("ic_method_for_flavors",
+ &(E->trace.ic_method_for_flavors),"0,0,nomax",m);
+ switch(E->trace.ic_method_for_flavors){
+ case 0: /* layer */
+ input_double("z_interface",&(E->trace.z_interface),"0.7",m);
+ break;
+ case 1: /* from grid in top n materials */
+ input_string("ictracer_grd_file",E->trace.ggrd_file,"",m); /* file from which to read */
+ input_int("ictracer_grd_layers",&(E->trace.ggrd_layers),"2",m); /* which top layers to use */
+ break;
+ default:
+ fprintf(stderr,"ic_method_for_flavors %i undefined\n",E->trace.ic_method_for_flavors);
+ parallel_process_termination();
+ break;
+ }
+
@@ -119,6 +136,7 @@
composition_input(E);
+
}
return;
@@ -162,7 +180,7 @@
void tracer_advection(struct All_variables *E)
{
-
+ if(E->control.verbose)
fprintf(E->trace.fpt,"STEP %d\n",E->monitor.solution_cycles);
/* advect tracers */
@@ -197,14 +215,15 @@
double convection_time,tracer_time;
double trace_fraction,total_time;
+ if(E->control.verbose){
+ fprintf(E->trace.fpt,"Number of times for all element search %d\n",E->trace.istat1);
- fprintf(E->trace.fpt,"Number of times for all element search %d\n",E->trace.istat1);
+ fprintf(E->trace.fpt,"Number of tracers sent to other processors: %d\n",E->trace.istat_isend);
- fprintf(E->trace.fpt,"Number of tracers sent to other processors: %d\n",E->trace.istat_isend);
+ fprintf(E->trace.fpt,"Number of times element columns are checked: %d \n",E->trace.istat_elements_checked);
+ }
- fprintf(E->trace.fpt,"Number of times element columns are checked: %d \n",E->trace.istat_elements_checked);
-
/* reset statistical counters */
E->trace.istat_isend=0;
@@ -215,6 +234,7 @@
/* compositional and error fraction data files */
//TODO: move
if (E->composition.on) {
+ if(E->control.verbose)
fprintf(E->trace.fpt,"Empty elements filled with old compositional "
"values: %d (%f percent)\n", E->trace.istat_iempty,
(100.0*E->trace.istat_iempty)/E->lmesh.nel);
@@ -529,9 +549,9 @@
reduce_tracer_arrays(E);
time_stat2=CPU_time0();
+ if(E->control.verbose)
+ fprintf(E->trace.fpt,"AA: time for find tracers: %f\n", time_stat2-time_stat1);
- fprintf(E->trace.fpt,"AA: time for find tracers: %f\n", time_stat2-time_stat1);
-
return;
}
@@ -1099,21 +1119,39 @@
int j, kk, number_of_tracers;
double rad;
+ switch(E->trace.ic_method_for_flavors){
+ case 0:
+ /* ic_method_for_flavors == 0 (layered structure) */
+ /* any tracer above z_interface is of flavor 0 */
+ /* any tracer below z_interface is of flavor 1 */
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+
+ number_of_tracers = E->trace.ntracers[j];
+ for (kk=1;kk<=number_of_tracers;kk++) {
+ rad = E->trace.basicq[j][2][kk];
+
+ if (rad<=E->trace.z_interface) E->trace.extraq[j][0][kk]=1.0;
+ if (rad>E->trace.z_interface) E->trace.extraq[j][0][kk]=0.0;
+ }
+ }
+ break;
- /* ic_method_for_flavors == 0 (layered structure) */
- /* any tracer above z_interface is of flavor 0 */
- /* any tracer below z_interface is of flavor 1 */
- if (E->trace.ic_method_for_flavors == 0) {
- for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ case 1: /* from grd in top n layers */
+#ifndef USE_GGRD
+ fprintf(stderr,"ic_method_for_flavors %i requires the ggrd routines from hc, -DUSE_GGRD\n",
+ E->trace.ic_method_for_flavors);
+ parallel_process_termination();
+#else
+ ggrd_init_tracer_flavors(E);
+#endif
+ break;
- number_of_tracers = E->trace.ntracers[j];
- for (kk=1;kk<=number_of_tracers;kk++) {
- rad = E->trace.basicq[j][2][kk];
- if (rad<=E->trace.z_interface) E->trace.extraq[j][0][kk]=1.0;
- if (rad>E->trace.z_interface) E->trace.extraq[j][0][kk]=0.0;
- }
- }
+ default:
+
+ fprintf(stderr,"ic_method_for_flavors %i undefined\n",E->trace.ic_method_for_flavors);
+ parallel_process_termination();
+ break;
}
return;
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-08 22:55:39 UTC (rev 7788)
@@ -167,11 +167,8 @@
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.PDEPV) /* "plasticity" */
+ visc_from_P(E,evisc);
if(E->viscosity.channel || E->viscosity.wedge)
@@ -269,7 +266,7 @@
imark = 0;
switch (E->viscosity.RHEOL) {
- case 1:
+ case 1: /* eta = N_0 exp( E * (T_0 - T)) */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nel;i++) {
l = E->mat[m][i];
@@ -296,7 +293,7 @@
}
break;
- case 2:
+ case 2: /* eta = N_0 exp(-T/T_0) */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nel;i++) {
l = E->mat[m][i];
@@ -323,37 +320,33 @@
}
break;
- case 3:
+ case 3: /* eta = N_0 exp(E/(T+T_0) - E/(1+T_0)) */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nel;i++) {
l = E->mat[m][i];
- tempa = E->viscosity.N0[l-1];
+ if(E->control.mat_control) /* switch moved up here TWB */
+ tempa = E->viscosity.N0[l-1] * E->VIP[m][i];
+ else
+ tempa = E->viscosity.N0[l-1];
j = 0;
for(kk=1;kk<=ends;kk++) {
- TT[kk] = E->T[m][E->ien[m][i].node[kk]];
- zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
+ TT[kk] = E->T[m][E->ien[m][i].node[kk]];
}
for(jj=1;jj<=vpts;jj++) {
temp=0.0;
- zzz=0.0;
- for(kk=1;kk<=ends;kk++) {
- TT[kk]=max(TT[kk],zero);
- temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
- zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+ for(kk=1;kk<=ends;kk++) { /* took out
+ computation of
+ depth, not
+ needed TWB */
+ TT[kk]=max(TT[kk],zero);
+ temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
}
-
- if(E->control.mat_control==0)
- EEta[m][ (i-1)*vpts + jj ] = tempa*
- exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
- - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
-
- if(E->control.mat_control==1)
- EEta[m][ (i-1)*vpts + jj ] = tempa*E->VIP[m][i]*
- exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
- - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
+ EEta[m][ (i-1)*vpts + jj ] = tempa*
+ exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
+ - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
}
}
break;
@@ -363,7 +356,11 @@
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nel;i++) {
l = E->mat[m][i];
- tempa = E->viscosity.N0[l-1];
+ if(E->control.mat_control) /* moved this up here TWB */
+ tempa = E->viscosity.N0[l-1] * E->VIP[m][i];
+ else
+ tempa = E->viscosity.N0[l-1];
+
j = 0;
for(kk=1;kk<=ends;kk++) {
@@ -389,25 +386,16 @@
so: E->viscosity.E = Ea/R*DT ; E->viscosity.Z = Va/R*DT
p = zzz and E->viscosity.T = T0/DT */
-
- if(E->control.mat_control==0)
- EEta[m][ (i-1)*vpts + jj ] = tempa*
- exp( (E->viscosity.E[l-1] + E->viscosity.Z[l-1]*zzz )
- / (E->viscosity.T[l-1]+temp) );
-
-
-
- if(E->control.mat_control==1)
- EEta[m][ (i-1)*vpts + jj ] = tempa*E->VIP[m][i]*
- exp( (E->viscosity.E[l-1] + E->viscosity.Z[l-1]*zzz )
- / (E->viscosity.T[l-1]+temp) );
-
+ EEta[m][ (i-1)*vpts + jj ] = tempa*
+ exp( (E->viscosity.E[l-1] + E->viscosity.Z[l-1]*zzz )
+ / (E->viscosity.T[l-1]+temp) );
+
}
}
break;
- case 5:
+ case 5: /* this still needs to be documented, who wrote this? */
/* same as rheol 3, except alternative margin, VIP, formulation */
for(m=1;m<=E->sphere.caps_per_proc;m++)
@@ -458,8 +446,41 @@
break;
+ case 6: /* eta = N_0 exp(E(T_0-T) + (1-z) Z_0 ) */
+
+ for(m=1;m <= E->sphere.caps_per_proc;m++)
+ for(i=1;i <= nel;i++) {
+ l = E->mat[m][i];
+ if(E->control.mat_control)
+ tempa = E->viscosity.N0[l-1] * E->VIP[m][i];
+ else
+ tempa = E->viscosity.N0[l-1];
+ j = 0;
+
+ for(kk=1;kk<=ends;kk++) {
+ TT[kk] = E->T[m][E->ien[m][i].node[kk]];
+ zz[kk] = (1.0 - E->sx[m][3][E->ien[m][i].node[kk]]);
+ }
+
+ for(jj=1;jj <= vpts;jj++) {
+ temp=0.0;zzz=0.0;
+ for(kk=1;kk <= ends;kk++) {
+ TT[kk]=max(TT[kk],zero);
+ temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
+ zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+ }
+ //fprintf(stderr,"N0 %11g T %11g T0 %11g E %11g 1-z %11g Z %11g\n",
+ // tempa,temp,E->viscosity.T[l-1],E->viscosity.E[l-1], zzz ,E->viscosity.Z[l-1] );
+ EEta[m][ (i-1)*vpts + jj ] = tempa*
+ exp( E->viscosity.E[l-1]*(E->viscosity.T[l-1] - temp) +
+ zzz * E->viscosity.Z[l-1] );
+ }
+ }
+ break;
+
+
}
return;
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-08 22:55:39 UTC (rev 7788)
@@ -659,11 +659,9 @@
int comp_nd; /* whether to output composition at nodes */
-#ifdef USE_GZDIR
- /* flags only used by GZDIR */
+ /* flags used by GZDIR */
int gzdir_vtkio,gzdir_vtkbase_init,gzdir_vtkbase_save;
float *gzdir_vtkbase;
-#endif
};
Modified: mc/3D/CitcomS/trunk/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-08-08 16:41:20 UTC (rev 7787)
+++ mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-08-08 22:55:39 UTC (rev 7788)
@@ -27,6 +27,7 @@
*/
+
/* forward declaration */
struct All_variables;
@@ -67,7 +68,12 @@
int ic_method_for_flavors;
double z_interface;
+
+ char ggrd_file[255]; /* for grd input */
+ int ggrd_layers;
+
+
/* statistical parameters */
int istat_ichoice[13][5];
int istat_isend;
More information about the cig-commits
mailing list