[cig-commits] r7814 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Mon Aug 13 20:33:22 PDT 2007
Author: becker
Date: 2007-08-13 20:33:21 -0700 (Mon, 13 Aug 2007)
New Revision: 7814
Modified:
mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
mc/3D/CitcomS/trunk/lib/Full_version_dependent.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/Pan_problem_misc_functions.c
mc/3D/CitcomS/trunk/lib/Problem_related.c
mc/3D/CitcomS/trunk/lib/Process_buoyancy.c
mc/3D/CitcomS/trunk/lib/Regional_sphere_related.c
mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
mc/3D/CitcomS/trunk/lib/Tracer_setup.c
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
- added coor=2 option (for full citcom) to assign finer radial node
spacing to top and lower layers of shell. The
coor_refine=0.1,0.15,0.1,0.2 parameters specify the radius fraction
of the bottom layer [0], the fraction of the nodes in this layer
[1], the top layer fraction [2], and the top layer node fraction
[3]. I.e. the defaults will put 15% of all nz nodes into the 10%
lower layer, 20% in the top 10% upper layer, and the rest in
between.
- renamed gzipped output version with sub-directory storage ascii-gz
- built in restart facilities for temperature and tracers when using
ascii-gz I/O with vtkio != 2
- added a composition 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
- added tracer_enriched option for internal heating. If tracer = on
and tracer_enriched = on, will reader Q0_enriched and vary the element heat production
between Q0 for C = 0 and Q0_enriched for C = 1. I.e. this only works
if C varies between 0 and 1.
- added an option to write from all processros to a single VTK file,
if ascii-gz is activated, and vtkio = 2. The VTK output is of the
"legacy", serial, single-file type, and requires that all processors see the same
filesystem.
This will lead to a bottleneck for large # of CPU computations as
each processor has to wait til the previous is done.
More efficient I/O should be possible by using the distributed
storage version of VTK, but I have no clue how this works. Anyone?
Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-08-13 22:55:41 UTC (rev 7813)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-08-14 03:33:21 UTC (rev 7814)
@@ -505,8 +505,19 @@
Q += Q0.Q[i] * exp(-Q0.lambda[i] * (E->monitor.elapsed_time+Q0.t_offset));
*/
+ /* heat production */
Q = E->control.Q0;
+ /* should we add a compositional contribution? */
+ if(E->control.tracer_enriched){
+
+ /* Q = Q0 for C = 0, Q = Q0ER for C = 1, and linearly in
+ between */
+ Q *= (1.0 - E->composition.comp_el[m][el]);
+ Q += E->composition.comp_el[m][el] * E->control.Q0ER;
+ }
+
+
/* construct residual from this information */
Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-08-13 22:55:41 UTC (rev 7813)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-08-14 03:33:21 UTC (rev 7814)
@@ -30,6 +30,8 @@
#include "global_defs.h"
#include "parallel_related.h"
+void get_r_spacing_fine(double *, double , double ,int ,double , double ,double , double,struct All_variables *);
+
/* Setup global mesh parameters */
void full_global_derived_values(E)
struct All_variables *E;
@@ -147,28 +149,39 @@
rr = (double *) malloc((E->mesh.noz+1)*sizeof(double));
RR = (double *) malloc((E->mesh.noz+1)*sizeof(double));
- if(E->control.coor==1) {
- sprintf(output_file,"%s",E->control.coor_file);
- fp1=fopen(output_file,"r");
- if (fp1 == NULL) {
- fprintf(E->fp,"(Nodal_mesh.c #1) Cannot open %s\n",output_file);
- exit(8);
- }
- fscanf(fp1,"%s %d",a,&i);
- for (k=1;k<=E->mesh.noz;k++) {
- fscanf(fp1,"%d %f",&nn,&tt1);
- rr[k]=tt1;
- }
- fclose(fp1);
- }
- else {
+ switch(E->control.coor){
+ case 2:
+ /* higher radial spacing in top and bottom fractions */
+ get_r_spacing_fine(rr, (double)E->sphere.ri,(double)E->sphere.ro,
+ E->mesh.noz,(double)E->control.coor_refine[0] ,
+ (double)E->control.coor_refine[1] ,
+ (double)E->control.coor_refine[2] ,
+ (double)E->control.coor_refine[3], E);
+ break;
+ case 1: /* read nodal radii from file */
+ sprintf(output_file,"%s",E->control.coor_file);
+ fp1=fopen(output_file,"r");
+ if (fp1 == NULL) {
+ fprintf(E->fp,"(Nodal_mesh.c #1) Cannot open %s\n",output_file);
+ exit(8);
+ }
+ fscanf(fp1,"%s %d",a,&i);
+ for (k=1;k<=E->mesh.noz;k++) {
+ fscanf(fp1,"%d %f",&nn,&tt1);
+ rr[k]=tt1;
+ }
+
+ fclose(fp1);
+ break;
+ default:
/* generate uniform mesh in radial direction */
dr = (E->sphere.ro-E->sphere.ri)/(E->mesh.noz-1);
-
+
for (k=1;k<=E->mesh.noz;k++) {
rr[k] = E->sphere.ri + (k-1)*dr;
}
+ break;
}
for (i=1;i<=E->lmesh.noz;i++) {
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-13 22:55:41 UTC (rev 7813)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-14 03:33:21 UTC (rev 7814)
@@ -291,6 +291,13 @@
input_int("levels",&(E->mesh.levels),"0",m);
input_int("coor",&(E->control.coor),"0",m);
+ if(E->control.coor == 2){ /* refinement */
+ E->control.coor_refine[0] = 0.10; /* bottom 10% */
+ E->control.coor_refine[1] = 0.15; /* get 15% of the nodes */
+ E->control.coor_refine[2] = 0.10; /* top 10% */
+ E->control.coor_refine[3] = 0.20; /* get 20% of the nodes */
+ input_float_vector("coor_refine",4,E->control.coor_refine,m);
+ }
input_string("coor_file",E->control.coor_file,"",m);
input_int("nprocx",&(E->parallel.nprocx),"1",m);
@@ -374,6 +381,8 @@
/* data section */
input_float("Q0",&(E->control.Q0),"0.0",m);
+ /* Q0_enriched gets read in Tracer_setup.c */
+
input_float("gravacc",&(E->data.grav_acc),"9.81",m);
input_float("thermexp",&(E->data.therm_exp),"3.0e-5",m);
input_float("cp",&(E->data.Cp),"1200.0",m);
@@ -1223,6 +1232,7 @@
void output_finalize(struct All_variables *E)
{
+ char message[255];
if (E->fp)
fclose(E->fp);
@@ -1232,6 +1242,20 @@
if (E->fp_out)
fclose(E->fp_out);
+ /*
+ remove VTK geo file in case we used that for IO (we'll only do
+ this for one processor, since this IO option requires shared
+ filesystems anyway
+ */
+ if((E->parallel.me == 0) && (E->output.gzdir.vtk_io == 2)&&
+ (strcmp(E->output.format, "ascii-gz") == 0)){
+ /* delete the vtk geo pre-file */
+ snprintf(message,255,"rm -f %s/vtk_geo",
+ E->control.data_dir);
+ system(message);
+ /* close the log */
+ fclose(E->output.gzdir.vtk_fp);
+ }
}
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2007-08-13 22:55:41 UTC (rev 7813)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2007-08-14 03:33:21 UTC (rev 7814)
@@ -49,6 +49,7 @@
void output_horiz_avg(struct All_variables *, int);
void output_tracer(struct All_variables *, int);
void output_pressure(struct All_variables *, int);
+FILE* output_open_mode(char *, char *);
extern void parallel_process_termination();
extern void heat_flux(struct All_variables *);
@@ -66,11 +67,16 @@
/* gzdir type of I/O */
if(strcmp(E->output.format, "ascii-gz") == 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!) */
+ /*
+ vtk_io = 1: write files for post-processing into VTK
+ vtk_io = 2: write legacy VTK file straight
+
+ */
+ input_int("gzdir_vtkio",&(E->output.gzdir.vtk_io),"0",m);
+ 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",
- // E->output.gzdir_vtkio,E->output.gzdir_vtkbase_save);
+ // E->output.gzdir.vtk_io,E->output.gzdir.vtk_base_save);
}
}
@@ -119,13 +125,18 @@
FILE* output_open(char *filename)
{
+ return output_open_mode(filename,"w");
+}
+FILE* output_open_mode(char *filename, char *mode)
+{
FILE *fp1;
/* if filename is empty, output to stderr. */
if (*filename) {
- fp1 = fopen(filename,"w");
+ fp1 = fopen(filename,mode);
if (!fp1) {
- fprintf(stderr,"Cannot open file '%s'\n",filename);
+ fprintf(stderr,"Cannot open file '%s' for '%s'\n",
+ filename,mode);
parallel_process_termination();
}
}
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-08-13 22:55:41 UTC (rev 7813)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-08-14 03:33:21 UTC (rev 7814)
@@ -30,31 +30,61 @@
/*
-this version uses gzipped, ascii output to subdirectories
+this version uses gzipped, ascii output to subdirectories for the
+ascii-gz option
+if, additionally, gzdir.vtk_io = 1, will write different format files
+ for post-processing into VTK
+
+ gzdir.vtk_io = 2, will try to write VTK straight (experimental)
+
+ the VTK output is the "legacy" type, requires that
+ all processors see the same filesystem, and will
+ likely lead to a bottleneck for larg CPU
+ computations as each processor has to wait til the
+ previous is done.
+
TWB
*/
-
#ifdef USE_GZDIR
+
+//#define ASCII_DEBUG
+
#include <zlib.h>
-
+#define BE_WERROR {myerror(E,"write error be output");}
+#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "element_definitions.h"
#include "global_defs.h"
#include "parsing.h"
#include "output.h"
+/* Big endian crap */
+#include <string.h>
+#include <malloc.h>
+void be_flipit(void *, void *, size_t );
+void be_flip_byte_order(void *, size_t );
+int be_is_little_endian(void);
+int be_write_float_to_file(float *, int , FILE *);
+int be_write_int_to_file(int *, int , FILE *);
+void myfprintf(FILE *,char *);
+
+/* */
+FILE* output_open_mode(char *, char *);
+void get_vtk_filename(char *,int,struct All_variables *);
+
+
gzFile *gzdir_output_open(char *,char *);
void gzdir_output(struct All_variables *, int );
void gzdir_output_comp_nd(struct All_variables *, int);
void gzdir_output_comp_el(struct All_variables *, int);
void gzdir_output_coord(struct All_variables *);
void gzdir_output_mat(struct All_variables *);
-void gzdir_output_velo(struct All_variables *, int);
+void gzdir_output_velo_temp(struct All_variables *, int);
void gzdir_output_visc_prepare(struct All_variables *, float **);
void gzdir_output_visc(struct All_variables *, int);
void gzdir_output_surf_botm(struct All_variables *, int);
@@ -87,32 +117,30 @@
void gzdir_output(struct All_variables *E, int cycles)
{
char output_dir[255];
-
if (cycles == 0) {
/* initial I/O */
gzdir_output_coord(E);
/*gzdir_output_mat(E);*/
}
+
/*
make a new directory for all the other output
- I thought I could just have the first proc do this and sync, but
- the syncing didn't work ?!
-
+ (all procs need to do that, because we might be using a local tmp
+ dir)
+
*/
- //if(E->parallel.me == 0){
- /* make a directory */
- snprintf(output_dir,255,"%s/%d",
- E->control.data_dir,cycles);
- mkdatadir(output_dir);
- // }
- /* and wait for the other jobs */
- // parallel_process_sync();
+ /* make a directory */
+ snprintf(output_dir,255,"%s/%d",E->control.data_dir,cycles);
+ mkdatadir(output_dir);
+
/* output */
- gzdir_output_velo(E, cycles);
+ gzdir_output_velo_temp(E, cycles); /* don't move this around,
+ else new VTK output won't
+ work */
gzdir_output_visc(E, cycles);
gzdir_output_surf_botm(E, cycles);
@@ -132,7 +160,8 @@
gzdir_output_horiz_avg(E, cycles);
if(E->control.tracer){
- if(E->output.tracer || (cycles == E->advection.max_timesteps))
+ if(E->output.tracer ||
+ (cycles == E->advection.max_timesteps))
gzdir_output_tracer(E, cycles);
}
@@ -166,208 +195,451 @@
void gzdir_output_coord(struct All_variables *E)
{
- int i, j, offset;
- char output_file[255];
- float locx[3];
- gzFile *fp1;
- /*
- don't use data file name
- */
- snprintf(output_file,255,"%s/coord.%d.gz",
- E->control.data_dir,E->parallel.me);
- fp1 = gzdir_output_open(output_file,"w");
+ int i, j, offset,ix[9],out;
+ char output_file[255],ostring[255],message[255];
+ float x[3];
+ gzFile *gz1;
+ FILE *fp1;
+
+ /* for dealing with several processors */
+ MPI_Status mpi_stat;
+ int mpi_rc, mpi_inmsg, mpi_success_message = 1;
+ if(E->output.gzdir.vtk_io == 2){
- /* nodal coordinates */
- for(j=1;j<=E->sphere.caps_per_proc;j++) {
- gzprintf(fp1,"%3d %7d\n",j,E->lmesh.nno);
- for(i=1;i<=E->lmesh.nno;i++)
- gzprintf(fp1,"%.6e %.6e %.6e\n",
- E->sx[j][1][i],E->sx[j][2][i],E->sx[j][3][i]);
- }
+ parallel_process_sync(E);
+ /*
+ new VTK output, write directly to a VTK file
+ this will get deleted from
- gzclose(fp1);
+ output_finalize(struct All_variables *E)
-
- if(E->output.gzdir_vtkio){
- /*
- output of Cartesian coordinates and element connectivitiy for
- vtk visualization
-
*/
- /*
- nodal coordinates in Cartesian
- */
- snprintf(output_file,255,"%s/vtk_ecor.%d.gz",
- E->control.data_dir,E->parallel.me);
- fp1 = gzdir_output_open(output_file,"w");
+ E->output.gzdir.vtk_ocount = -1;
+ snprintf(output_file,255,"%s/vtk_geo",E->control.data_dir);
+ if(E->parallel.me == 0){
+
+ /* log times */
+ snprintf(message,255,
+ "%s/vtk_time.log",E->control.data_dir);
+ E->output.gzdir.vtk_fp = output_open_mode(message,"w");
+
+
+ /* start geo file */
+ fp1 = output_open_mode(output_file,"w");
+ myfprintf(fp1,"# vtk DataFile Version 2.0\n");
+ myfprintf(fp1,"model name, extra info\n");
+#ifdef ASCII_DEBUG
+ myfprintf(fp1,"ASCII\n");
+#else
+ myfprintf(fp1,"BINARY\n");
+#endif
+ myfprintf(fp1,"DATASET UNSTRUCTURED_GRID\n");
+ sprintf(message,"POINTS %i float\n", /* total number of nodes */
+ E->lmesh.nno * E->parallel.nproc *
+ E->sphere.caps_per_proc);
+ myfprintf(fp1,message);
+ }else{
+ /* if not first, wait for previous */
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), 0, MPI_COMM_WORLD, &mpi_stat);
+ /* open for append */
+ fp1 = output_open_mode(output_file,"a");
+ }
+ out = 0;
+ /* write nodal coordinate to file, big endian */
for(j=1;j <= E->sphere.caps_per_proc;j++) {
for(i=1;i <= E->lmesh.nno;i++) {
- rtp2xyz(E->sx[j][3][i],E->sx[j][1][i],E->sx[j][2][i],locx);
- gzprintf(fp1,"%9.6f %9.6f %9.6f\n",
- locx[0],locx[1],locx[2]);
+ x[0]=E->x[j][1][i];x[1]=E->x[j][2][i];x[2]=E->x[j][3][i];
+ if(be_write_float_to_file(x,3,fp1) != 3)BE_WERROR;
+ out++;
}
}
- gzclose(fp1);
+ fclose(fp1);fflush(fp1); /* close file and flush buffer */
+ if(E->parallel.me < E->parallel.nproc-1){/* send to next if not last*/
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 0, MPI_COMM_WORLD);
+ }
/*
- connectivity for all elements
+ node numbers for all the elements
*/
+ parallel_process_sync(E);
+ if(E->parallel.me == 0){
+ fp1 = output_open_mode(output_file,"a");
+ j = E->parallel.nproc * E->lmesh.nel *
+ E->sphere.caps_per_proc; /* total number of elements */
+ sprintf(message,"CELLS %i %i\n", /* number of elements
+ total number of int entries
+
+ */
+ j,j*(enodes[E->mesh.nsd]+1));
+ myfprintf(fp1,message);
+ }else{
+ /* if not first, wait for previous */
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), 0, MPI_COMM_WORLD, &mpi_stat);
+ fp1 = output_open_mode(output_file,"a");
+ }
+ /*
+ write CELL element nodes
+ */
+ if(enodes[E->mesh.nsd] != 8)
+ myerror(E,"vtk error, only eight node hexes supported");
offset = E->lmesh.nno * E->parallel.me - 1;
- snprintf(output_file,255,"%s/vtk_econ.%d.gz",
- E->control.data_dir,E->parallel.me);
- fp1 = gzdir_output_open(output_file,"w");
+ ix[0] = enodes[E->mesh.nsd];
for(j=1;j <= E->sphere.caps_per_proc;j++) {
for(i=1;i <= E->lmesh.nel;i++) {
- gzprintf(fp1,"%2i\t",enodes[E->mesh.nsd]);
- if(enodes[E->mesh.nsd] != 8){
- gzprintf(stderr,"gzdir: Output: error, only eight node hexes supported");
- parallel_process_termination();
- }
/*
need to add offset according to the processor for global
node numbers
*/
- gzprintf(fp1,"%6i %6i %6i %6i %6i %6i %6i %6i\n",
- E->ien[j][i].node[1]+offset,E->ien[j][i].node[2]+offset,
- E->ien[j][i].node[3]+offset,E->ien[j][i].node[4]+offset,
- E->ien[j][i].node[5]+offset,E->ien[j][i].node[6]+offset,
- E->ien[j][i].node[7]+offset,E->ien[j][i].node[8]+offset);
+ ix[1]= E->ien[j][i].node[1]+offset;ix[2] = E->ien[j][i].node[2]+offset;
+ ix[3]= E->ien[j][i].node[3]+offset;ix[4] = E->ien[j][i].node[4]+offset;
+ ix[5]= E->ien[j][i].node[5]+offset;ix[6] = E->ien[j][i].node[6]+offset;
+ ix[7]= E->ien[j][i].node[7]+offset;ix[8] = E->ien[j][i].node[8]+offset;
+ if(be_write_int_to_file(ix,9,fp1)!=9)BE_WERROR;
}
}
- gzclose(fp1);
- } /* end vtkio */
+ fclose(fp1);fflush(fp1); /* close file and flush buffer */
+ if(E->parallel.me < E->parallel.nproc-1)
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 0, MPI_COMM_WORLD);
+ /*
+ add the cell type flags
+ */
+ parallel_process_sync(E);
+ if(E->parallel.me == 0){
+ fp1 = output_open_mode(output_file,"a");
+ j=E->parallel.nproc*E->lmesh.nel*E->sphere.caps_per_proc;
+ sprintf(message,"CELL_TYPES %i\n",j); /* number of elements*/
+ myfprintf(fp1,message);
+ ix[0] = 12;
+ for(i=0;i<j;i++)
+ if(be_write_int_to_file(ix,1,fp1)!=1)BE_WERROR;
+ fclose(fp1);fflush(fp1); /* close file and flush buffer */
+ report(E,"vtk_io: vtk geometry done\n");
+ }
+ /* done straight VTK */
+ }else{
+ /*
+
+ either zipped regular, or old VTK type for post-processing
+
+ */
+ /*
+ don't use data file name
+ */
+ snprintf(output_file,255,"%s/coord.%d.gz",
+ E->control.data_dir,E->parallel.me);
+ gz1 = gzdir_output_open(output_file,"w");
+
+ /* nodal coordinates */
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ gzprintf(gz1,"%3d %7d\n",j,E->lmesh.nno);
+ for(i=1;i<=E->lmesh.nno;i++)
+ gzprintf(gz1,"%.6e %.6e %.6e\n",
+ E->sx[j][1][i],E->sx[j][2][i],E->sx[j][3][i]);
+ }
+
+ gzclose(gz1);
+ if(E->output.gzdir.vtk_io){
+ /*
+
+ output of Cartesian coordinates and element connectivitiy for
+ vtk visualization
+
+ */
+ /*
+ nodal coordinates in Cartesian
+ */
+ snprintf(output_file,255,"%s/vtk_ecor.%d.gz",
+ E->control.data_dir,E->parallel.me);
+ gz1 = gzdir_output_open(output_file,"w");
+ for(j=1;j <= E->sphere.caps_per_proc;j++) {
+ for(i=1;i <= E->lmesh.nno;i++) {
+ gzprintf(gz1,"%9.6f %9.6f %9.6f\n", /* cartesian nodal coordinates */
+ E->x[j][1][i],E->x[j][2][i],E->x[j][3][i]);
+ }
+ }
+ gzclose(gz1);
+ /*
+ connectivity for all elements
+ */
+ offset = E->lmesh.nno * E->parallel.me - 1;
+ snprintf(output_file,255,"%s/vtk_econ.%d.gz",
+ E->control.data_dir,E->parallel.me);
+ gz1 = gzdir_output_open(output_file,"w");
+ for(j=1;j <= E->sphere.caps_per_proc;j++) {
+ for(i=1;i <= E->lmesh.nel;i++) {
+ gzprintf(gz1,"%2i\t",enodes[E->mesh.nsd]);
+ if(enodes[E->mesh.nsd] != 8){
+ gzprintf(stderr,"gzdir: Output: error, only eight node hexes supported");
+ parallel_process_termination();
+ }
+ /*
+ need to add offset according to the processor for global
+ node numbers
+ */
+ gzprintf(gz1,"%6i %6i %6i %6i %6i %6i %6i %6i\n",
+ E->ien[j][i].node[1]+offset,E->ien[j][i].node[2]+offset,
+ E->ien[j][i].node[3]+offset,E->ien[j][i].node[4]+offset,
+ E->ien[j][i].node[5]+offset,E->ien[j][i].node[6]+offset,
+ E->ien[j][i].node[7]+offset,E->ien[j][i].node[8]+offset);
+ }
+ }
+ gzclose(gz1);
+ } /* end vtkio = 1 */
+ }
-
return;
}
+/*
-void gzdir_output_visc(struct All_variables *E, int cycles)
-{
- int i, j;
- char output_file[255];
- gzFile *fp1;
- int lev = E->mesh.levmax;
+this needs to be called first before any of the other stuff if VTK
+straight output is chosen
- snprintf(output_file,255,
- "%s/%d/visc.%d.%d.gz", E->control.data_dir,
- cycles,E->parallel.me, cycles);
- fp1 = gzdir_output_open(output_file,"w");
- for(j=1;j<=E->sphere.caps_per_proc;j++) {
- gzprintf(fp1,"%3d %7d\n",j,E->lmesh.nno);
- for(i=1;i<=E->lmesh.nno;i++)
- gzprintf(fp1,"%.4e\n",E->VI[lev][j][i]);
- }
-
- gzclose(fp1);
- return;
-}
-
-
-void gzdir_output_velo(struct All_variables *E, int cycles)
+*/
+void gzdir_output_velo_temp(struct All_variables *E, int cycles)
{
int i, j, k,os;
- char output_file[255],output_file2[255];
+ char output_file[255],output_file2[255],message[255];
float cvec[3];
gzFile *gzout;
- /*
-
+ FILE *fp1;
+ /* for dealing with several processors */
+ MPI_Status mpi_stat;
+ int mpi_rc;
+ int mpi_inmsg, mpi_success_message = 1;
- temperatures are printed along with velocities for old type of
- output
- if VTK is selected, will generate a separate temperature file
-
- */
- if(E->output.gzdir_vtkio) {
- /*
- for VTK, only print temperature
- */
- snprintf(output_file2,255,"%s/%d/t.%d.%d",
- E->control.data_dir,
- cycles,E->parallel.me,cycles);
- }else{ /* vel + T */
- snprintf(output_file2,255,"%s/%d/velo.%d.%d",
- E->control.data_dir,cycles,
- E->parallel.me,cycles);
- }
- snprintf(output_file,255,"%s.gz",output_file2); /* add the .gz */
+ if(E->output.gzdir.vtk_io){ /* all VTK modes need basis vectors */
- gzout = gzdir_output_open(output_file,"w");
- gzprintf(gzout,"%d %d %.5e\n",
- cycles,E->lmesh.nno,E->monitor.elapsed_time);
- for(j=1; j<= E->sphere.caps_per_proc;j++) {
- gzprintf(gzout,"%3d %7d\n",j,E->lmesh.nno);
- if(E->output.gzdir_vtkio){
- /* VTK */
- for(i=1;i<=E->lmesh.nno;i++)
- gzprintf(gzout,"%.6e\n",E->T[j][i]);
- } else {
- /* old */
- for(i=1;i<=E->lmesh.nno;i++)
- gzprintf(gzout,"%.6e %.6e %.6e %.6e\n",
- E->sphere.cap[j].V[1][i],
- E->sphere.cap[j].V[2][i],
- E->sphere.cap[j].V[3][i],E->T[j][i]);
+ os = E->lmesh.nno*9;
+
+ if((!E->output.gzdir.vtk_base_init) ||(!E->output.gzdir.vtk_base_save)){
+ if(!E->output.gzdir.vtk_base_init)
+ E->output.gzdir.vtk_base = (float *)safe_malloc(sizeof(float)*os*E->sphere.caps_per_proc);
+ for(k=0,j=1;j <= E->sphere.caps_per_proc;j++,k += os) {
+ for(i=1;i <= E->lmesh.nno;i++,k += 9){
+ /* cartesian basis vectors at theta, phi */
+ calc_cbase_at_tp(E->sx[j][1][i],E->sx[j][2][i],(E->output.gzdir.vtk_base+k));
+ }
+ }
+ E->output.gzdir.vtk_base_init = 1;
}
}
- gzclose(gzout);
- if(E->output.gzdir_vtkio){
+
+ if(E->output.gzdir.vtk_io == 2){
+ parallel_process_sync(E);
+
/*
- cartesian velocity output
+
+ straight VTK output, always need to start here!
+
*/
- os = E->lmesh.nno*9;
+ E->output.gzdir.vtk_ocount++;
+ get_vtk_filename(output_file,E->output.gzdir.vtk_ocount,E);
/*
- get base vectors if first pass or if we're not saving
+
+ start with temperature
+
*/
- if((!E->output.gzdir_vtkbase_init) ||
- (!E->output.gzdir_vtkbase_save)){
- if(!E->output.gzdir_vtkbase_init){
- /* allocate */
- E->output.gzdir_vtkbase = (float *)
- safe_malloc(sizeof(float)*os*E->sphere.caps_per_proc);
+ if(E->parallel.me == 0){
+ /* copy geo file over to start out vtk file */
+ snprintf(output_file2,255,"cp %s/vtk_geo %s",
+ E->control.data_dir,output_file);
+ system(output_file2);
+ /* should we do something to check if this has worked? */
+
+ /* write a time log */
+ fprintf(E->output.gzdir.vtk_fp,"%12i %12i %12.6e %s\n",
+ E->output.gzdir.vtk_ocount,cycles,E->monitor.elapsed_time,output_file);
+
+ /* */
+ fp1 = output_open_mode(output_file,"a");
+ sprintf(message,"POINT_DATA %i\n",E->lmesh.nno*E->parallel.nproc*E->sphere.caps_per_proc);
+ myfprintf(fp1,message);
+ myfprintf(fp1,"SCALARS temperature float 1\n");
+ myfprintf(fp1,"LOOKUP_TABLE default\n");
+ }else{
+ /* if not first, wait for previous */
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), 7, MPI_COMM_WORLD, &mpi_stat);
+ /* open for append */
+ fp1 = output_open_mode(output_file,"a");
+ }
+ for(j=1; j<= E->sphere.caps_per_proc;j++) /* print the temperatures */
+ for(i=1;i<=E->lmesh.nno;i++){
+ cvec[0] = E->T[j][i];
+ if(be_write_float_to_file(cvec,1,fp1)!=1)BE_WERROR;
}
- for(k=0,j=1;j <= E->sphere.caps_per_proc;j++,k += os) {
- for(i=1;i <= E->lmesh.nno;i++,k += 9){
- /* cartesian basis vectors at theta, phi */
- calc_cbase_at_tp(E->sx[j][1][i],
- E->sx[j][2][i],
- (E->output.gzdir_vtkbase+k));
- }
- }
- E->output.gzdir_vtkbase_init = 1;
+ fclose(fp1);fflush(fp1); /* close file and flush buffer */
+ if(E->parallel.me < E->parallel.nproc-1){
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 7, MPI_COMM_WORLD);
+ }else{
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, 0, 6, MPI_COMM_WORLD); /* tell m=0 to go ahead */
}
/*
- write Cartesian velocities to file
+ velocities second
*/
- snprintf(output_file,255,"%s/%d/vtk_v.%d.%d.gz",
- E->control.data_dir,cycles,E->parallel.me,cycles);
- gzout = gzdir_output_open(output_file,"w");
+ if(E->parallel.me == 0){
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, E->parallel.nproc-1 , 6, MPI_COMM_WORLD, &mpi_stat);
+ fp1 = output_open_mode(output_file,"a"); /* append velocities */
+ sprintf(message,"VECTORS velocity float\n");myfprintf(fp1,message);
+ }else{
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), 5, MPI_COMM_WORLD, &mpi_stat);
+ fp1 = output_open_mode(output_file,"a");
+ }
for(k=0,j=1;j <= E->sphere.caps_per_proc;j++,k += os) {
for(i=1;i<=E->lmesh.nno;i++,k += 9) {
- /* convert r,theta,phi vector to x,y,z at base location */
- convert_pvec_to_cvec(E->sphere.cap[j].V[3][i],
- E->sphere.cap[j].V[1][i],
- E->sphere.cap[j].V[2][i],
- (E->output.gzdir_vtkbase+k),cvec);
- /* output of cartesian vector */
- gzprintf(gzout,"%10.4e %10.4e %10.4e\n",
- cvec[0],cvec[1],cvec[2]);
+ convert_pvec_to_cvec(E->sphere.cap[j].V[3][i],E->sphere.cap[j].V[1][i],E->sphere.cap[j].V[2][i],
+ (E->output.gzdir.vtk_base+k),cvec);
+ if(be_write_float_to_file(cvec,3,fp1)!=3)BE_WERROR;
}
}
+ fclose(fp1);fflush(fp1); /* close file and flush buffer */
+ if(E->parallel.me < E->parallel.nproc-1){
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 5, MPI_COMM_WORLD);
+ }else{
+ fprintf(stderr,"vtk_io: geo, temp, & vel writtend to %s\n",output_file);
+ }
+ /* new VTK velo and temp done */
+ }else{
+ /*
+
+ modified zipped output
+
+ */
+ /*
+
+
+ temperatures are printed along with velocities for old type of
+ output
+
+ if VTK is selected, will generate a separate temperature file
+
+ */
+ if(E->output.gzdir.vtk_io) {
+ /*
+ for VTK, only print temperature
+ */
+ snprintf(output_file2,255,"%s/%d/t.%d.%d",
+ E->control.data_dir,
+ cycles,E->parallel.me,cycles);
+ }else{ /* vel + T */
+ snprintf(output_file2,255,"%s/%d/velo.%d.%d",
+ E->control.data_dir,cycles,
+ E->parallel.me,cycles);
+ }
+ snprintf(output_file,255,"%s.gz",output_file2); /* add the .gz */
+
+ gzout = gzdir_output_open(output_file,"w");
+ gzprintf(gzout,"%d %d %.5e\n",
+ cycles,E->lmesh.nno,E->monitor.elapsed_time);
+ for(j=1; j<= E->sphere.caps_per_proc;j++) {
+ gzprintf(gzout,"%3d %7d\n",j,E->lmesh.nno);
+ if(E->output.gzdir.vtk_io){
+ /* VTK */
+ for(i=1;i<=E->lmesh.nno;i++)
+ gzprintf(gzout,"%.6e\n",E->T[j][i]);
+ } else {
+ /* old */
+ for(i=1;i<=E->lmesh.nno;i++)
+ gzprintf(gzout,"%.6e %.6e %.6e %.6e\n",
+ E->sphere.cap[j].V[1][i],
+ E->sphere.cap[j].V[2][i],
+ E->sphere.cap[j].V[3][i],E->T[j][i]);
+ }
+ }
gzclose(gzout);
-
+ if(E->output.gzdir.vtk_io){
+ /*
+ write Cartesian velocities to file
+ */
+ snprintf(output_file,255,"%s/%d/vtk_v.%d.%d.gz",
+ E->control.data_dir,cycles,E->parallel.me,cycles);
+ gzout = gzdir_output_open(output_file,"w");
+ for(k=0,j=1;j <= E->sphere.caps_per_proc;j++,k += os) {
+ for(i=1;i<=E->lmesh.nno;i++,k += 9) {
+ /* convert r,theta,phi vector to x,y,z at base location */
+ convert_pvec_to_cvec(E->sphere.cap[j].V[3][i],
+ E->sphere.cap[j].V[1][i],
+ E->sphere.cap[j].V[2][i],
+ (E->output.gzdir.vtk_base+k),cvec);
+ /* output of cartesian vector */
+ gzprintf(gzout,"%10.4e %10.4e %10.4e\n",
+ cvec[0],cvec[1],cvec[2]);
+ }
+ }
+ gzclose(gzout);
+
+ }
+ } /* end gzipped and old VTK out */
+ if(E->output.gzdir.vtk_io){ /* all VTK modes */
/* free memory */
- if(!E->output.gzdir_vtkbase_save)
- free(E->output.gzdir_vtkbase);
+ if(!E->output.gzdir.vtk_base_save)
+ free(E->output.gzdir.vtk_base);
}
+ return;
+}
+/*
+ viscosity
+*/
+void gzdir_output_visc(struct All_variables *E, int cycles)
+{
+ int i, j;
+ char output_file[255];
+ gzFile *gz1;
+ FILE *fp1;
+ int lev = E->mesh.levmax;
+ float ftmp;
+ /* for dealing with several processors */
+ MPI_Status mpi_stat;
+ int mpi_rc;
+ int mpi_inmsg, mpi_success_message = 1;
+
+ if(E->output.gzdir.vtk_io != 2){ /* old output */
+ snprintf(output_file,255,
+ "%s/%d/visc.%d.%d.gz", E->control.data_dir,
+ cycles,E->parallel.me, cycles);
+ gz1 = gzdir_output_open(output_file,"w");
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ gzprintf(gz1,"%3d %7d\n",j,E->lmesh.nno);
+ for(i=1;i<=E->lmesh.nno;i++)
+ gzprintf(gz1,"%.4e\n",E->VI[lev][j][i]);
+ }
+
+ gzclose(gz1);
+ }else{
+ parallel_process_sync(E);
+
+ /* new legacy VTK */
+ get_vtk_filename(output_file,E->output.gzdir.vtk_ocount,E);
+ if(E->parallel.me == 0){
+ fp1 = output_open_mode(output_file,"a");
+ myfprintf(fp1,"SCALARS log10(visc) float 1\n");
+ myfprintf(fp1,"LOOKUP_TABLE default\n");
+ }else{
+ /* if not first, wait for previous */
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), 0, MPI_COMM_WORLD, &mpi_stat);
+ /* open for append */
+ fp1 = output_open_mode(output_file,"a");
+ }
+ for(j=1; j<= E->sphere.caps_per_proc;j++)
+ for(i=1;i<=E->lmesh.nno;i++){
+ ftmp = log10(E->VI[lev][j][i]);
+ if(fabs(ftmp) < 5e-7)ftmp = 0.0;
+ if(be_write_float_to_file(&ftmp,1,fp1)!=1)BE_WERROR;
+ }
+ fclose(fp1);fflush(fp1); /* close file and flush buffer */
+ if(E->parallel.me < E->parallel.nproc-1){
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 0, MPI_COMM_WORLD);
+ }
+ }
return;
}
+
void gzdir_output_surf_botm(struct All_variables *E, int cycles)
{
int i, j, s;
@@ -543,24 +815,47 @@
void gzdir_output_pressure(struct All_variables *E, int cycles)
{
int i, j;
+ float ftmp;
char output_file[255];
- gzFile *fp1;
+ gzFile *gz1;
+ FILE *fp1;
+ /* for dealing with several processors */
+ MPI_Status mpi_stat;
+ int mpi_rc;
+ int mpi_inmsg, mpi_success_message = 1;
- 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");
-
- 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);
- for(i=1;i<=E->lmesh.nno;i++)
- gzprintf(fp1,"%.6e\n",E->NP[j][i]);
+ if(E->output.gzdir.vtk_io != 2){ /* old */
+ snprintf(output_file,255,"%s/%d/pressure.%d.%d.gz", E->control.data_dir,cycles,
+ E->parallel.me, cycles);
+ gz1 = gzdir_output_open(output_file,"w");
+ gzprintf(gz1,"%d %d %.5e\n",cycles,E->lmesh.nno,E->monitor.elapsed_time);
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ gzprintf(gz1,"%3d %7d\n",j,E->lmesh.nno);
+ for(i=1;i<=E->lmesh.nno;i++)
+ gzprintf(gz1,"%.6e\n",E->NP[j][i]);
+ }
+ gzclose(gz1);
+ }else{/* new legacy VTK */
+ parallel_process_sync(E);
+ get_vtk_filename(output_file,E->output.gzdir.vtk_ocount,E);
+ if(E->parallel.me == 0){
+ fp1 = output_open_mode(output_file,"a");
+ myfprintf(fp1,"SCALARS pressure float 1\n");
+ myfprintf(fp1,"LOOKUP_TABLE default\n");
+ }else{
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), 0, MPI_COMM_WORLD, &mpi_stat);
+ fp1 = output_open_mode(output_file,"a");
+ }
+ for(j=1; j<= E->sphere.caps_per_proc;j++)
+ for(i=1;i<=E->lmesh.nno;i++){
+ ftmp = E->NP[j][i];
+ if(be_write_float_to_file(&ftmp,1,fp1)!=1)BE_WERROR;
+ }
+ fclose(fp1);fflush(fp1); /* close file and flush buffer */
+ if(E->parallel.me < E->parallel.nproc-1){
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 0, MPI_COMM_WORLD);
+ }
}
-
- gzclose(fp1);
-
return;
}
@@ -606,32 +901,54 @@
void gzdir_output_comp_nd(struct All_variables *E, int cycles)
{
- int i, j;
- char output_file[255];
- gzFile *fp1;
-
+ int i, j;
+ char output_file[255];
+ gzFile *gz1;
+ FILE *fp1;
+ float ftmp;
+ /* for dealing with several processors */
+ MPI_Status mpi_stat;
+ int mpi_rc;
+ int mpi_inmsg, mpi_success_message = 1;
+
+ if(E->output.gzdir.vtk_io != 2){
snprintf(output_file,255,"%s/%d/comp_nd.%d.%d.gz",
- E->control.data_dir,
- cycles,
- E->parallel.me, cycles);
- fp1 = gzdir_output_open(output_file,"w");
-
+ E->control.data_dir,cycles,
+ E->parallel.me, cycles);
+ gz1 = gzdir_output_open(output_file,"w");
for(j=1;j<=E->sphere.caps_per_proc;j++) {
- gzprintf(fp1,"%3d %7d %.5e %.5e %.5e\n",
- j, E->lmesh.nel,
- E->monitor.elapsed_time,
- E->composition.initial_bulk_composition,
- E->composition.bulk_composition);
-
- for(i=1;i<=E->lmesh.nno;i++) {
- gzprintf(fp1,"%.6e\n",
- E->composition.comp_node[j][i]);
- }
-
+ gzprintf(gz1,"%3d %7d %.5e %.5e %.5e\n",
+ j, E->lmesh.nel,
+ E->monitor.elapsed_time,
+ E->composition.initial_bulk_composition,
+ E->composition.bulk_composition);
+ for(i=1;i<=E->lmesh.nno;i++) {
+ gzprintf(gz1,"%.6e\n",E->composition.comp_node[j][i]);
+ }
}
-
- gzclose(fp1);
- return;
+ gzclose(gz1);
+ }else{/* new legacy VTK */
+ parallel_process_sync(E);
+ get_vtk_filename(output_file,E->output.gzdir.vtk_ocount,E);
+ if(E->parallel.me == 0){
+ fp1 = output_open_mode(output_file,"a");
+ myfprintf(fp1,"SCALARS composition float 1\n");
+ myfprintf(fp1,"LOOKUP_TABLE default\n");
+ }else{
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), 0, MPI_COMM_WORLD, &mpi_stat);
+ fp1 = output_open_mode(output_file,"a");
+ }
+ for(j=1; j<= E->sphere.caps_per_proc;j++)
+ for(i=1;i<=E->lmesh.nno;i++){
+ ftmp = E->composition.comp_node[j][i];
+ if(be_write_float_to_file(&ftmp,1,fp1)!=1)BE_WERROR;
+ }
+ fclose(fp1);fflush(fp1); /* close file and flush buffer */
+ if(E->parallel.me < E->parallel.nproc-1){
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT, (E->parallel.me+1), 0, MPI_COMM_WORLD);
+ }
+ }
+ return;
}
@@ -678,7 +995,7 @@
ii = E->monitor.solution_cycles_init;
- if(E->output.gzdir_vtkio) { /* VTK I/O */
+ if(E->output.gzdir.vtk_io) { /* VTK I/O */
snprintf(output_file,255,"%s/%d/t.%d.%d",
E->control.data_dir_old,
ii,E->parallel.me,ii);
@@ -696,7 +1013,7 @@
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 */
+ if(E->output.gzdir.vtk_io) { /* 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");
@@ -784,4 +1101,135 @@
}
+
+
+void get_vtk_filename(char *output_file,int cycles,struct All_variables *E)
+{
+ snprintf(output_file,255,"%s/d.%08i.vtk",E->control.data_dir,cycles);
+}
+
+
+
+
+/*
+
+
+big endian I/O (needed for vtk)
+
+
+*/
+
+/*
+
+write the x[n] array to file, making sure it is written big endian
+
+*/
+int be_write_float_to_file(float *x, int n, FILE *out)
+{
+ int i,nout;
+ static size_t len = sizeof(float);
+ size_t bsize;
+ float ftmp;
+#ifdef ASCII_DEBUG
+ for(i=0;i<n;i++)
+ fprintf(out,"%11g ",x[i]);
+ fprintf(out,"\n");
+ nout = n;
+#else
+ /*
+ do we need to flip?
+ */
+ if(be_is_little_endian()){
+ nout = 0;
+ for(i=0;i < n;i++){
+ ftmp = x[i];
+ be_flip_byte_order((void *)(&ftmp),len);
+ nout += fwrite(&ftmp,len,(size_t)1,out); /* write to file */
+ }
+ }else{ /* operate on x */
+ nout = fwrite(x,len,(size_t)n,out); /* write to file */
+ }
+#endif
+ return nout;
+}
+int be_write_int_to_file(int *x, int n, FILE *out)
+{
+ int i,nout;
+ static size_t len = sizeof(int);
+ size_t bsize;
+ int itmp;
+#ifdef ASCII_DEBUG
+ for(i=0;i<n;i++)
+ fprintf(out,"%11i ",x[i]);
+ fprintf(out,"\n");
+ nout = n;
+#else
+ /*
+ do we need to flip?
+ */
+ if(be_is_little_endian()){
+ nout = 0;
+ for(i=0;i < n;i++){
+ itmp = x[i];
+ be_flip_byte_order((void *)(&itmp),len);
+ nout += fwrite(&itmp,len,(size_t)1,out); /* write to file */
+ }
+ }else{ /* operate on x */
+ nout = fwrite(x,len,(size_t)n,out); /* write to file */
+ }
+#endif
+ return nout;
+}
+
+
+/* does this make a difference? nope, didn't, and why would it */
+void myfprintf(FILE *out,char *string)
+{
+#ifdef ASCII_DEBUG
+ fprintf(out,string);
+#else
+ fwrite(string, sizeof(char), strlen(string), out);
+#endif
+}
+
+int be_is_little_endian(void)
+{
+ static const unsigned long a = 1;
+ return *(const unsigned char *)&a;
+}
+
+/*
+
+
+flip endian-ness
+
+
+*/
+/*
+
+flip endianness of x
+
+*/
+void be_flip_byte_order(void *x, size_t len)
+{
+ void *copy;
+ int i;
+ copy = (void *)malloc(len); /* don't check here for speed */
+ memcpy(copy,x,len);
+ be_flipit(x,copy,len);
+ free(copy);
+}
+
+/* this should not be called with (i,i,size i) */
+void be_flipit(void *d, void *s, size_t len)
+{
+ unsigned char *dest = d;
+ unsigned char *src = s;
+ src += len - 1;
+ for (; len; len--)
+ *dest++ = *src--;
+}
+
+
+#undef BE_WERROR
#endif /* gzdir switch */
Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-13 22:55:41 UTC (rev 7813)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-14 03:33:21 UTC (rev 7814)
@@ -56,6 +56,8 @@
void convert_pvec_to_cvec(float ,float , float , float *,float *);
void *safe_malloc (size_t );
void myerror(struct All_variables *,char *);
+void get_r_spacing_fine(double *, double , double ,int ,double , double ,double , double,
+ struct All_variables *);
@@ -473,3 +475,49 @@
E->parallel.me,message);
parallel_process_termination();
}
+
+/*
+
+
+
+attempt to space rr[1...nz] such that bfrac*nz nodes will be within the lower
+brange fraction of (ro-ri), and similar for the top layer
+
+*/
+void get_r_spacing_fine(double *rr, double ri, double ro,
+ int nz,double brange, double bfrac,
+ double trange, double tfrac,
+ struct All_variables *E)
+{
+ int k,klim,nb,nt,nm;
+ double drb,dr0,drt,dr,drm,range,r,mrange;
+ range = ro-ri; /* original range */
+ mrange = 1 - brange - trange;
+ if(mrange <= 0)
+ myerror(E,"get_r_spacing_fine: bottom and top range too large");
+
+ brange *= range; /* bottom */
+ trange *= range; /* top */
+ mrange *= range; /* middle */
+
+ nb = nz * bfrac;
+ nt = nz * tfrac;
+ nm = nz-nb-nt;
+ if((nm < 1)||(nt < 2)||(nb < 2))
+ myerror(E,"get_r_spacing_fine: refinement out of bounds");
+
+ drb = brange/(nb-1);
+ drt = trange/(nt-1);
+ drm = mrange / (nm + 1);
+
+ for(r=ri,k=1;k<=nb;k++,r+=drb){
+ rr[k] = r;
+ }
+ klim = nz-nt+1;
+ for(r=r-drb+drm;k < klim;k++,r+=drm){
+ rr[k] = r;
+ }
+ for(;k <= nz;k++,r+=drt){
+ rr[k] = r;
+ }
+}
Modified: mc/3D/CitcomS/trunk/lib/Problem_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Problem_related.c 2007-08-13 22:55:41 UTC (rev 7813)
+++ mc/3D/CitcomS/trunk/lib/Problem_related.c 2007-08-14 03:33:21 UTC (rev 7814)
@@ -74,7 +74,7 @@
#ifdef USE_GZDIR /* gzdir output */
if(strcmp(E->output.format, "ascii-gz") == 0){
- if(E->output.gzdir_vtkio)
+ if(E->output.gzdir.vtk_io)
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
Modified: mc/3D/CitcomS/trunk/lib/Process_buoyancy.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Process_buoyancy.c 2007-08-13 22:55:41 UTC (rev 7813)
+++ mc/3D/CitcomS/trunk/lib/Process_buoyancy.c 2007-08-14 03:33:21 UTC (rev 7814)
@@ -35,9 +35,9 @@
#include "element_definitions.h"
#include "global_defs.h"
+#include <math.h> /* for sqrt */
-
void post_processing(struct All_variables *E)
{
return;
Modified: mc/3D/CitcomS/trunk/lib/Regional_sphere_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_sphere_related.c 2007-08-13 22:55:41 UTC (rev 7813)
+++ mc/3D/CitcomS/trunk/lib/Regional_sphere_related.c 2007-08-14 03:33:21 UTC (rev 7814)
@@ -161,7 +161,7 @@
free ((void *)fi1[lev] );
}
- }
+}
else if(E->control.coor==0) {
@@ -243,6 +243,9 @@
free ((void *)SX[1]);
free ((void *)tt);
free ((void *)ff);
+}else{
+ myerror(E,"regional_sphere_related: coor setting not implemented");
+
}
return;
Modified: mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2007-08-13 22:55:41 UTC (rev 7813)
+++ mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2007-08-14 03:33:21 UTC (rev 7814)
@@ -29,6 +29,7 @@
#include "global_defs.h"
#include "parallel_related.h"
+void get_r_spacing_fine(double *, double , double ,int ,double , double ,double , double,struct All_variables *);
/* Setup global mesh parameters */
void regional_global_derived_values(E)
@@ -161,44 +162,52 @@
noz=E->mesh.noz;
- if(E->control.coor==1) {
- sprintf(output_file,"%s",E->control.coor_file);
- fp1=fopen(output_file,"r");
- if (fp1 == NULL) {
- fprintf(E->fp,"(Nodal_mesh.c #1) Cannot open %s\n",output_file);
- exit(8);
- }
-
- fscanf(fp1,"%s %d",a,&i);
- for(i=1;i<=nox;i++)
+ if(E->control.coor==1) { /* get nodal levels from file */
+ sprintf(output_file,"%s",E->control.coor_file);
+ fp1=fopen(output_file,"r");
+ if (fp1 == NULL) {
+ fprintf(E->fp,"(Nodal_mesh.c #1) Cannot open %s\n",output_file);
+ exit(8);
+ }
+
+ fscanf(fp1,"%s %d",a,&i);
+ for(i=1;i<=nox;i++)
fscanf(fp1,"%d %f",&nn,&tt1);
-
- fscanf(fp1,"%s %d",a,&i);
- for(i=1;i<=noy;i++)
+
+ fscanf(fp1,"%s %d",a,&i);
+ for(i=1;i<=noy;i++)
fscanf(fp1,"%d %f",&nn,&tt1);
-
- fscanf(fp1,"%s %d",a,&i);
- for (k=1;k<=E->mesh.noz;k++) {
+
+ fscanf(fp1,"%s %d",a,&i);
+ for (k=1;k<=E->mesh.noz;k++) {
fscanf(fp1,"%d %f",&nn,&tt1);
rr[k]=tt1;
- }
- E->sphere.ri = rr[1];
- E->sphere.ro = rr[E->mesh.noz];
-
- fclose(fp1);
-
- }
-
- else {
- dr = (E->sphere.ro-E->sphere.ri)/(E->mesh.noz-1);
- for (k=1;k<=E->mesh.noz;k++) {
+ }
+ E->sphere.ri = rr[1];
+ E->sphere.ro = rr[E->mesh.noz];
+
+ fclose(fp1);
+
+ } else if(E->control.coor==0) {
+ /* default: regular node spacing */
+ dr = (E->sphere.ro-E->sphere.ri)/(E->mesh.noz-1);
+ for (k=1;k<=E->mesh.noz;k++) {
rr[k] = E->sphere.ri + (k-1)*dr;
- }
-
}
+ } else if(E->control.coor==2){
+ /* higher radial spacing in top and bottom fractions */
+ get_r_spacing_fine(rr, (double)E->sphere.ri,(double)E->sphere.ro,
+ E->mesh.noz,(double)E->control.coor_refine[0] ,
+ (double)E->control.coor_refine[1] ,
+ (double)E->control.coor_refine[2] ,
+ (double)E->control.coor_refine[3],E);
+ } else {
+ myerror(E,"regional_version_dependent: coor mode not implemented");
+ }
+
for (i=1;i<=E->lmesh.noz;i++) {
k = E->lmesh.nzs+i-1;
RR[i] = rr[k];
Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-08-13 22:55:41 UTC (rev 7813)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-08-14 03:33:21 UTC (rev 7814)
@@ -82,9 +82,22 @@
void tracer_input(struct All_variables *E)
{
void full_tracer_input();
+ char message[100];
int m=E->parallel.me;
- input_boolean("tracer",&(E->control.tracer),"0",m);
+ input_boolean("tracer",&(E->control.tracer),"off",m);
+ input_boolean("tracer_enriched",
+ &(E->control.tracer_enriched),"off",m);
+ if(E->control.tracer_enriched){
+ if(!E->control.tracer) /* check here so that we can get away
+ with only one if statement in
+ Advection_diffusion */
+ myerror(E,"need to switch on tracers for tracer_enriched");
+ input_float("Q0_enriched",&(E->control.Q0ER),"0.0",m);
+ snprintf(message,100,"using compositionally enriched heating: C = 0: %g C = 1: %g",
+ E->control.Q0,E->control.Q0ER);
+ report(E,message);
+ }
if(E->control.tracer) {
/* tracer_ic_method=0 (random generated array) */
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-13 22:55:41 UTC (rev 7813)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-14 03:33:21 UTC (rev 7814)
@@ -496,6 +496,7 @@
int pseudo_free_surf;
int tracer;
+ int tracer_enriched;
double theta_min, theta_max, fi_min, fi_max;
float start_age;
@@ -524,6 +525,7 @@
float VBYbotval;
int coor;
+ float coor_refine[4];
char coor_file[100];
int lith_age;
@@ -541,6 +543,8 @@
float TBCbotval;
float Q0;
+ float Q0ER;
+
float jrelax;
int precondition;
@@ -625,6 +629,15 @@
float timedir;
};
+ /* for gzdir */
+struct gzd_struc{
+ int vtk_io,vtk_base_init,vtk_base_save,
+ vtk_ocount;
+ float *vtk_base;
+ FILE *vtk_fp;
+
+};
+
struct Output {
char format[20]; /* ascii or hdf5 */
char optional[1000]; /* comma-delimited list of objects to output */
@@ -660,8 +673,7 @@
/* flags used by GZDIR */
- int gzdir_vtkio,gzdir_vtkbase_init,gzdir_vtkbase_save;
- float *gzdir_vtkbase;
+ struct gzd_struc gzdir;
};
More information about the cig-commits
mailing list