[cig-commits] r4124 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Mon Jul 31 14:02:15 PDT 2006
Author: becker
Date: 2006-07-31 14:02:14 -0700 (Mon, 31 Jul 2006)
New Revision: 4124
Modified:
mc/1D/hc/trunk/Makefile
mc/1D/hc/trunk/README.TXT
mc/1D/hc/trunk/ggrd_grdtrack_util.c
mc/1D/hc/trunk/ggrd_grdtrack_util.h
mc/1D/hc/trunk/ggrd_velinterpol.c
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_extract_sh_layer.c
mc/1D/hc/trunk/hc_filenames.h
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/hc_output.c
mc/1D/hc/trunk/hc_polsol.c
mc/1D/hc/trunk/hc_solve.c
mc/1D/hc/trunk/main.c
Log:
Added traction functionality.
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/Makefile 2006-07-31 21:02:14 UTC (rev 4124)
@@ -51,7 +51,7 @@
PREM_DEFINES = -DPREM_MODEL_FILE=\"$(PWD)/prem/prem.dat\"
PREM_INCS = prem.h
#
-# GMT grd handling, need to define USE_GMT4 for version 4.0
+# GMT grd handling
#
GGRD_SRCS = ggrd_velinterpol.c ggrd_readgrds.c ggrd_grdtrack_util.c
GGRD_OBJS = $(ODIR)/ggrd_velinterpol.o $(ODIR)/ggrd_readgrds.o $(ODIR)/ggrd_grdtrack_util.o
Modified: mc/1D/hc/trunk/README.TXT
===================================================================
--- mc/1D/hc/trunk/README.TXT 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/README.TXT 2006-07-31 21:02:14 UTC (rev 4124)
@@ -1,3 +1,54 @@
+INSTALLATION
+
+The code should compile on any basic UNIX/Linux system. See Makefile
+for comments.
+
+
+HC CODE usage
+
+Described in the help page that is displayed for "hc -h" as below. See
+also the benchmark scripts in cig_bench.tgz. This tar file will expand
+into a cig subdirectory particularly "run_benchmark" there shows how
+to run all CIG benchmark tests.
+
+
+>>>
+
+
+hc - perform Hager & O'Connell flow computation
+
+This code can compute velocities, tractions, and geoid for simple
+density distributions and plate velocities using the semi-analytical
+approach of Hager & O'Connell (1981). This particular implementation
+illustrates one possible way to combine the HC solver routines.
+
+Based on code by Brad Hager, Richard O'Connell, and Bernhard Steinberger.
+
+density anomaly options:
+-dens name use name as a SH density anomaly model (dens.sh.dat)
+ This expects density anomalies in % of PREM in DT format.
+-ds additional density anomaly scaling factor (0.2)
+
+Earth model options:
+-prem name set Earth model to name (/home/becker/progs/src/hc-svn/prem/prem.dat)
+-vf name set viscosity structure filename to name (visc.dat)
+ This file is in non_dim_radius viscosity[Pas] format
+
+boundary condition options:
+-fs perform free slip computation, else no slip or plates (OFF)
+-ns perform no slip computation, else try to read plate velocities (OFF)
+-pvel name set surface velocity file to name (pvel.sh.dat)
+ This file is based on a DT expansion of cm/yr velocity fields.
+
+solution procedure and I/O options:
+-ng do not compute and print the geoid (OFF)
+-pptsol print pol[6] and tor[2] solution vectors (OFF)
+-px print the spatial solution to file (OFF)
+-str compute srr,srt,srp tractions [MPa] instead of velocities [cm/yr] (default: vel)
+-v -vv -vvv: verbosity levels (0)
+
+<<<
+
SPHERICAL HARMONICS FORMAT
All single layer spherical harmonics are in the fully normalized,
@@ -3,5 +54,5 @@
physical convention, e.g. Dahlen & Tromp (1998).
-
+
All spherical harmonics files have an SH_HEADER
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2006-07-31 21:02:14 UTC (rev 4124)
@@ -17,6 +17,7 @@
#define ONEEIGHTYOVERPI 57.295779513082320876798154814105
#endif
#include <math.h>
+#include <string.h>
/*
@@ -108,8 +109,8 @@
// if(change_z_sign) /* reverse logic */
//g->zlevels_are_negative = (g->zlevels_are_negative)?(FALSE):(TRUE);
if(verbose){
- fprintf(stderr,"ggrd_grdtrack_init_general: initialized from %s, %s.\n",
- grdfile,(is_three)?("3-D"):("1-D"));
+ fprintf(stderr,"ggrd_grdtrack_init_general: initialized from %s, %s, bcflag: %s.\n",
+ grdfile,(is_three)?("3-D"):("1-D"), gmt_edgeinfo_string);
if(is_three){
fprintf(stderr,"ggrd_grdtrack_init_general: depth file %s, %i levels, %s.\n",
depth_file,g->nz,
@@ -155,6 +156,9 @@
/*
+
+for 3-D spherical
+
interpolation wrapper, uses r, theta, phi input. return value and TRUE if success,
undefined and FALSE else
*/
@@ -195,6 +199,8 @@
return result;
}
/*
+ for 3-D cartesian
+
interpolation wrapper, uses x, y, z input. return value and TRUE if success,
undefined and FALSE else
*/
@@ -232,9 +238,13 @@
}
/*
- interpolation wrapper, uses theta, phi input.
- return value and TRUE if success,
- undefined and FALSE else
+
+for 2-D spherical
+
+interpolation wrapper, uses theta, phi input.
+return value and TRUE if success,
+undefined and FALSE else
+
*/
unsigned char ggrd_grdtrack_interpolate_tp(double t,double p,
struct ggrd_gt *g,
@@ -270,8 +280,44 @@
#endif
return result;
}
+
/*
+ for 2-D cartesian
+ interpolation wrapper, uses x,y input
+ return value and TRUE if success,
+ undefined and FALSE else
+*/
+unsigned char ggrd_grdtrack_interpolate_xy(double xin,double yin,
+ struct ggrd_gt *g,
+ double *value,
+ unsigned char verbose)
+{
+ double x[3];
+ unsigned char result;
+ if(!g->init){ /* this will not necessarily work */
+ fprintf(stderr,"ggrd_grdtrack_interpolate_xy: error, g structure not initialized\n");
+ return FALSE;
+ }
+ if(g->is_three){
+ fprintf(stderr,"ggrd_grdtrack_interpolate_xy: error, g structure is not 2-D\n");
+ return FALSE;
+ }
+ x[0] = xin;
+ x[1] = yin;
+ x[2] = 0.0;
+#ifdef USE_GMT4
+ result = ggrd_grdtrack_interpolate(x,FALSE,g->grd,g->f,g->edgeinfo,g->mm,g->z,g->nz,
+ value,verbose,&g->bcr);
+#else
+ result = ggrd_grdtrack_interpolate(x,FALSE,g->grd,g->f,g->edgeinfo,g->mm,g->z,g->nz,
+ value,verbose);
+#endif
+ return result;
+}
+
+/*
+
free structure
*/
@@ -382,18 +428,22 @@
init first edgeinfo
*/
-
GMT_boundcond_init (*edgeinfo);
- if ((edgeinfo_string+2)){
+ /* check if geographic */
+ if (strlen(edgeinfo_string)>2){
GMT_boundcond_parse (*edgeinfo, (edgeinfo_string+2));
- if ((*edgeinfo)[0].gn)
+ if ((*edgeinfo+0)->gn)
*geographic_in = TRUE;
else
*geographic_in = FALSE;
}else{
*geographic_in = FALSE;
}
-
+ if(verbose >= 2)
+ if(*geographic_in)
+ fprintf(stderr,"ggrd_grdtrack_init: detected geographic region from geostring: %s\n",
+ edgeinfo_string);
+
*z = (float *) GMT_memory
(VNULL, (size_t)1, sizeof(float), "ggrd_grdtrack_init");
@@ -586,8 +636,9 @@
return 10;
}
#endif
-
- /* prepare the boundaries */
+ /*
+ prepare the boundaries
+ */
GMT_boundcond_param_prep ((*grd+i), (*edgeinfo+i));
if(i == 0){
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.h
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.h 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.h 2006-07-31 21:02:14 UTC (rev 4124)
@@ -27,6 +27,14 @@
unsigned char ggrd_grdtrack_interpolate_xyz(double ,double ,double ,
struct ggrd_gt *,double *,
unsigned char);
+unsigned char ggrd_grdtrack_interpolate_xy(double ,double ,
+ struct ggrd_gt *,
+ double *,
+ unsigned char );
+unsigned char ggrd_grdtrack_interpolate_tp(double ,double ,
+ struct ggrd_gt *,
+ double *,
+ unsigned char );
void ggrd_grdtrack_free_gstruc(struct ggrd_gt *);
@@ -45,25 +53,8 @@
unsigned char ggrd_grdtrack_interpolate(double *, unsigned char , struct GRD_HEADER *, float *,
struct GMT_EDGEINFO *, int, float *, int , double *,unsigned char,
struct GMT_BCR *);
+int ggrd_grdtrack_init(double *, double *, double *, double *, float **, int *, char *, struct GRD_HEADER **, struct GMT_EDGEINFO **, char *, unsigned char *, int *, unsigned char, char *, float **, int *, unsigned char, unsigned char, unsigned char, struct GMT_BCR *);
-int ggrd_grdtrack_init(double *, double *,double *, double *, /* geographic bounds,
- set all to zero to
- get the whole range from the
- input grid files
- */
- float **, /* data */
- int *, /* size of data */
- char *, /* name, or prefix, of grd file with scalars */
- struct GRD_HEADER **,
- struct GMT_EDGEINFO **,
- char *,unsigned char *,
- int *, /* [4] array with padding (output) */
- unsigned char _d, char *, /* depth file name */
- float **, /* layers, pass as NULL */
- int *, /* number of layers */
- unsigned char , /* linear/cubic? */
- unsigned char ,unsigned char,
- struct GMT_BCR *);
#else
unsigned char ggrd_grdtrack_interpolate(double *, unsigned char , struct GRD_HEADER *, float *,
struct GMT_EDGEINFO *, int, float *, int , double *,unsigned char);
Modified: mc/1D/hc/trunk/ggrd_velinterpol.c
===================================================================
--- mc/1D/hc/trunk/ggrd_velinterpol.c 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/ggrd_velinterpol.c 2006-07-31 21:02:14 UTC (rev 4124)
@@ -28,7 +28,8 @@
#include "hc.h"
-int ggrd_find_vel_and_der(GGRD_CPREC *xloc,GGRD_CPREC time,GGRD_CPREC dtrange,
+int ggrd_find_vel_and_der(GGRD_CPREC *xloc,
+ GGRD_CPREC time,GGRD_CPREC dtrange,
struct ggrd_vel *v,int order,
hc_boolean icalc_der,
hc_boolean verbose,
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/hc.h 2006-07-31 21:02:14 UTC (rev 4124)
@@ -90,7 +90,14 @@
struct hc_sm{
HC_PREC u[6][4];
};
+/*
+
+parameter structure to allow for settings that are specific to the
+implementation and higher level than the calls to the subroutines
+
+
+ */
struct hc_parameters{
hc_boolean compressible; /* compressibility following Panasyuk
& Steinberger */
@@ -105,6 +112,7 @@
hc_boolean verbose; /* debugging output? (0,1,2,3,4...) */
hc_boolean sol_binary_out; /* binary or ASCII output of SH expansion */
hc_boolean print_spatial; /* print the spatial solution */
+ hc_boolean compute_geoid; /* compute and print the geoid */
int solution_mode; /* velocity or stress */
hc_boolean print_pt_sol; /* output of p[6] and t[2] vectors */
char visc_filename[HC_CHAR_LENGTH]; /* name of viscosity profile file */
@@ -225,6 +233,7 @@
velocity scale to change from input
[cm/yr] to internal
*/
+ HC_PREC stress_scale; /* to go to MPa */
HC_PREC r_cmb; /* radius of CMB */
/* Legendre functions */
double *plm;
@@ -242,7 +251,7 @@
#define HC_VEL 0 /* compute velocities
v_r v_t v_p
*/
-#define HC_STRESS -1 /* compute tractions
+#define HC_TRACTIONS 1 /* compute tractions
trr trt trp */
/*
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2006-07-31 21:02:14 UTC (rev 4124)
@@ -4,6 +4,7 @@
unsigned char ggrd_grdtrack_interpolate_rtp(double, double, double, struct ggrd_gt *, double *, unsigned char);
unsigned char ggrd_grdtrack_interpolate_xyz(double, double, double, struct ggrd_gt *, double *, unsigned char);
unsigned char ggrd_grdtrack_interpolate_tp(double, double, struct ggrd_gt *, double *, unsigned char);
+unsigned char ggrd_grdtrack_interpolate_xy(double, double, struct ggrd_gt *, double *, unsigned char);
void ggrd_grdtrack_free_gstruc(struct ggrd_gt *);
void ggrd_find_spherical_vel_from_rigid_cart_rot(double *, double *, double *, double *, double *);
int ggrd_grdtrack_init(double *, double *, double *, double *, float **, int *, char *, struct GRD_HEADER **, struct GMT_EDGEINFO **, char *, unsigned char *, int *, unsigned char, char *, float **, int *, unsigned char, unsigned char, unsigned char);
@@ -78,7 +79,7 @@
void hc_print_vector_label(double *, int, FILE *, char *);
void hc_print_matrix_label(double *, int, int, FILE *, char *);
void hc_print_vector_row(double *, int, FILE *);
-void hc_compute_solution_scaling_factors(struct hcs *, int, double *);
+void hc_compute_solution_scaling_factors(struct hcs *, int, double, double *);
void hc_print_poloidal_solution(struct sh_lms *, struct hcs *, int, char *, unsigned short, unsigned short);
void hc_print_toroidal_solution(double *, int, struct hcs *, int, char *, unsigned short);
/* hc_polsol.c */
Modified: mc/1D/hc/trunk/hc_extract_sh_layer.c
===================================================================
--- mc/1D/hc/trunk/hc_extract_sh_layer.c 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/hc_extract_sh_layer.c 2006-07-31 21:02:14 UTC (rev 4124)
@@ -44,15 +44,15 @@
default:
fprintf(stderr,"%s: usage\n%s sol.file layer [mode,%i] [short_format, %i]\n\n",
argv[0],argv[0],mode,short_format);
- fprintf(stderr,"extracts spherical harmonic solution from HC run\n");
+ fprintf(stderr,"extracts spherical harmonic solution x (vel or str) from HC run\n");
fprintf(stderr,"layer: 1...nset\n");
fprintf(stderr,"\tif ilayer=1..nset, will print one layer\n");
fprintf(stderr,"\tif ilayer=-1, will select nset\n");
fprintf(stderr,"\tif ilayer=-2, will print all layers\n");
fprintf(stderr,"mode: 1...3\n");
- fprintf(stderr,"\tif mode = 1, will print u_r \n");
- fprintf(stderr,"\tif mode = 2, will print u_pol u_tor \n");
- fprintf(stderr,"\tif mode = 3, will print u_r u_pol u_tor\n");
+ fprintf(stderr,"\tif mode = 1, will print x_r \n");
+ fprintf(stderr,"\tif mode = 2, will print x_pol x_tor \n");
+ fprintf(stderr,"\tif mode = 3, will print x_r x_pol x_tor\n");
exit(-1);
break;
@@ -101,21 +101,21 @@
case 1:
/* */
if(verbose)
- fprintf(stderr,"%s: printing u_r SHE at layer %i (depth: %g)\n",
+ fprintf(stderr,"%s: printing x_r SHE at layer %i (depth: %g)\n",
argv[0],ilayer,HC_Z_DEPTH(model->r[ilayer]));
sh_print_coefficients_to_file((sol+ilayer*3),shps,stdout,fac,FALSE,verbose);
break;
case 2:
/* */
if(verbose)
- fprintf(stderr,"%s: printing u_pol u_tor SHE at layer %i (depth: %g)\n",
+ fprintf(stderr,"%s: printing x_pol x_tor SHE at layer %i (depth: %g)\n",
argv[0],ilayer,HC_Z_DEPTH(model->r[ilayer]));
sh_print_coefficients_to_file((sol+ilayer*3+1),shps,stdout,fac,FALSE,verbose);
break;
case 3:
/* mode == 3 */
if(verbose)
- fprintf(stderr,"%s: printing u_r u_pol u_tor SHE at layer %i (depth: %g)\n",
+ fprintf(stderr,"%s: printing x_r x_pol x_tor SHE at layer %i (depth: %g)\n",
argv[0],ilayer,HC_Z_DEPTH(model->r[ilayer]));
sh_print_coefficients_to_file((sol+ilayer*3),shps,stdout,fac,FALSE,verbose);
break;
Modified: mc/1D/hc/trunk/hc_filenames.h
===================================================================
--- mc/1D/hc/trunk/hc_filenames.h 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/hc_filenames.h 2006-07-31 21:02:14 UTC (rev 4124)
@@ -13,12 +13,13 @@
#define HC_SOLOUT_FILE_ASCII "sol.dat" /* solution output files */
#define HC_SOLOUT_FILE_BINARY "sol.bin"
-#define HC_SPATIAL_SOLOUT_FILE "vsol" /* spatial solution
- output files, those will ahave
- stuff appended like 5.bin */
+#define HC_SPATIAL_SOLOUT_FILE "ssol" /* spatial solution output
+ files, those will ahave stuff
+ appended like 5.bin */
-#define HC_TORSOL_FILE "tsol.dat" /* toroidal solutions 1 and 2 as f(l,r) */
+#define HC_TORSOL_FILE "tsol.dat" /* toroidal solutions 1 and 2 as
+ f(l,r) */
#define HC_POLSOL_FILE "psol.dat"
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/hc_init.c 2006-07-31 21:02:14 UTC (rev 4124)
@@ -25,13 +25,13 @@
if false, plate velocities, else no
slip if free_slip is false as well
*/
+ p->compute_geoid = TRUE; /* compute the geoid? */
p->dens_anom_scale = 0.2; /* default density anomaly scaling to
go from PREM percent traveltime
anomalies to density anomalies */
p->verbose = 0; /* debugging output? (0,1,2,3,4...) */
p->sol_binary_out = TRUE; /* binary or ASCII output of SH expansion */
- p->solution_mode = HC_VEL; /* velocity, stress, or geoid */
-
+ p->solution_mode = HC_VEL; /* default: velocity */
p->print_pt_sol = FALSE;
p->print_spatial = FALSE; /* by default, only print the spectral solution */
/*
@@ -172,12 +172,15 @@
/*
constants
*/
- hc->timesc = HC_TIMESCALE_YR; /* timescale [yr]*/
+ hc->timesc = HC_TIMESCALE_YR; /* timescale [yr], like 1e6 yrs */
hc->visnor = 1e21; /* normalizing viscosity [Pas]*/
hc->gacc = 10.0e2; /* gravitational acceleration [cm/s2]*/
hc->g = 6.6742e-11; /* gravitational constant [Nm2/kg2]*/
+
/*
- radius of Earth in [m]
+
+ radius of Earth in [m]
+
*/
hc->re = hc->prem->r0;
if(fabs(hc->re - (HC_RE_KM * 1e3)) > 1e-7)
@@ -197,11 +200,26 @@
HC_ERROR("hc_init_constants","Earth model CMB radius appears off");
/*
+
+ derived quantities
+ */
+ /*
+
velocity scale if input is in [cm/yr], works out to be ~0.11
*/
hc->vel_scale = hc->re*PIOVERONEEIGHTY/hc->timesc*100;
+ /*
+
+ stress scaling, will later be divided by non-dim radius, to go
+ to MPa
+
+ */
+ hc->stress_scale = (PIOVERONEEIGHTY * hc->visnor / hc->secyr)/
+ (hc->timesc * 1e6);
+
+
init = TRUE;
}
@@ -223,29 +241,44 @@
*/
fprintf(stderr,"%s - perform Hager & O'Connell flow computation\n\n",
argv[0]);
- fprintf(stderr,"options:\n\n");
+ fprintf(stderr,"This code can compute velocities, tractions, and geoid for simple density distributions\n");
+ fprintf(stderr,"and plate velocities using the semi-analytical approach of Hager & O'Connell (1981).\n");
+ fprintf(stderr,"This particular implementation illustrates one possible way to combine the HC solver routines.\n\n");
+ fprintf(stderr,"Based on code by Brad Hager, Richard O'Connell, and Bernhard Steinberger.\n\n");
+
+ fprintf(stderr,"density anomaly options:\n");
fprintf(stderr,"-dens\tname\tuse name as a SH density anomaly model (%s)\n",
p->dens_filename);
- fprintf(stderr,"-ds\t\tdensity scaling (%g)\n",
+ fprintf(stderr,"\t\tThis expects density anomalies in %% of PREM in DT format.\n");
+ fprintf(stderr,"-ds\t\tadditional density anomaly scaling factor (%g)\n\n",
p->dens_anom_scale);
+
+ fprintf(stderr,"Earth model options:\n");
+ fprintf(stderr,"-prem\tname\tset Earth model to name (%s)\n",
+ p->prem_model_filename);
+ fprintf(stderr,"-vf\tname\tset viscosity structure filename to name (%s)\n",
+ p->visc_filename);
+ fprintf(stderr,"\t\tThis file is in non_dim_radius viscosity[Pas] format\n\n");
+
+ fprintf(stderr,"boundary condition options:\n");
fprintf(stderr,"-fs\t\tperform free slip computation, else no slip or plates (%s)\n",
hc_name_boolean(p->free_slip));
fprintf(stderr,"-ns\t\tperform no slip computation, else try to read plate velocities (%s)\n",
hc_name_boolean(p->vel_bc_zero));
- fprintf(stderr,"-prem\tname\tset Earth model to name (%s)\n",
- p->prem_model_filename);
- fprintf(stderr,"-pptsol\t\tprint pol[6] and tor[2] solution vectors (%s)\n",
- hc_name_boolean(p->print_pt_sol));
fprintf(stderr,"-pvel\tname\tset surface velocity file to name (%s)\n",
p->pvel_filename);
+ fprintf(stderr,"\t\tThis file is based on a DT expansion of cm/yr velocity fields.\n\n");
+
+ fprintf(stderr,"solution procedure and I/O options:\n");
+ fprintf(stderr,"-ng\t\tdo not compute and print the geoid (%s)\n",
+ hc_name_boolean(!p->compute_geoid));
+ fprintf(stderr,"-pptsol\t\tprint pol[6] and tor[2] solution vectors (%s)\n",
+ hc_name_boolean(p->print_pt_sol));
fprintf(stderr,"-px\t\tprint the spatial solution to file (%s)\n",
hc_name_boolean(p->print_spatial));
- fprintf(stderr,"-vf\tname\tset viscosity structure filename to name (%s)\n",
- p->visc_filename);
+ fprintf(stderr,"-str\t\tcompute srr,srt,srp tractions [MPa] instead of velocities [cm/yr] (default: vel)\n");
fprintf(stderr,"-v\t-vv\t-vvv: verbosity levels (%i)\n",
(int)(p->verbose));
- fprintf(stderr,"-visc\tname\tuse name as a viscosity structure (%s)\n",
- p->visc_filename);
fprintf(stderr,"\n\n");
exit(-1);
}else if(strcmp(argv[i],"-ds")==0){ /* density anomaly scaling factor */
@@ -262,6 +295,8 @@
solution
parameters */
hc_toggle_boolean(&p->print_pt_sol);
+ }else if(strcmp(argv[i],"-ng")==0){ /* do not compute geoid */
+ hc_toggle_boolean(&p->compute_geoid);
}else if(strcmp(argv[i],"-vf")==0){ /* viscosity filename */
hc_advance_argument(&i,argc,argv);
strncpy(p->visc_filename,argv[i],HC_CHAR_LENGTH);
@@ -277,6 +312,8 @@
}else if(strcmp(argv[i],"-visc")==0){ /* viscosity filename */
hc_advance_argument(&i,argc,argv);
strncpy(p->visc_filename,argv[i],HC_CHAR_LENGTH);
+ }else if(strcmp(argv[i],"-str")==0){ /* compute tractions */
+ p->solution_mode = HC_TRACTIONS;
}else if(strcmp(argv[i],"-v")==0){ /* verbosities */
p->verbose = 1;
}else if(strcmp(argv[i],"-vv")==0){ /* verbosities */
Modified: mc/1D/hc/trunk/hc_output.c
===================================================================
--- mc/1D/hc/trunk/hc_output.c 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/hc_output.c 2006-07-31 21:02:14 UTC (rev 4124)
@@ -29,14 +29,14 @@
number of solution sets of ntype solutions
*/
nradp2 = hc->nrad+2;
- /*
-
- scale to cm/yr, or other values for stress solutions
-
- */
- hc_compute_solution_scaling_factors(hc,sol_mode,fac);
for(i=os=0;i < nradp2;i++,os += ntype){
/*
+
+ scale to cm/yr, or MPa for stress solutions
+
+ */
+ hc_compute_solution_scaling_factors(hc,sol_mode,hc->r[i],fac);
+ /*
write parameters, convert radius to depth in [km]
*/
sh_print_parameters_to_file((sol+os),ntype,i,nradp2,
@@ -47,16 +47,24 @@
write the set of coefficients in D&T convention
*/
- sh_print_coefficients_to_file((sol+os),ntype,out,fac,binary,verbose);
+ sh_print_coefficients_to_file((sol+os),ntype,out,fac,
+ binary,verbose);
if(verbose >= 2){
switch(sol_mode){
case HC_VEL:
- fprintf(stderr,"hc_print_spectral_solution: z: %8.3f |r|: %11.3e |pol|: %11.3e |tor|: %11.3e (scale: %g cm/yr)\n",
+ fprintf(stderr,"hc_print_spectral_solution: z: %8.3f vel: |r|: %11.3e |pol|: %11.3e |tor|: %11.3e (scale: %g cm/yr)\n",
HC_Z_DEPTH(hc->r[i]),sqrt(sh_total_power((sol+os))),
sqrt(sh_total_power((sol+os+1))),
sqrt(sh_total_power((sol+os+2))),
fac[0]/11.1194926644559);
break;
+ case HC_TRACTIONS:
+ fprintf(stderr,"hc_print_spectral_solution: z: %8.3f str: |r|: %11.3e |pol|: %11.3e |tor|: %11.3e (scale: %g MPa)\n",
+ HC_Z_DEPTH(hc->r[i]),sqrt(sh_total_power((sol+os))),
+ sqrt(sh_total_power((sol+os+1))),
+ sqrt(sh_total_power((sol+os+2))),
+ fac[0]/(0.553073278428428/hc->r[i]));
+ break;
default:
fprintf(stderr,"hc_print_spectral_solution: sol mode %i undefined\n",sol_mode);
exit(-1);
@@ -96,7 +104,8 @@
will also write the corresponding depth layers to dfilename
*/
-void hc_print_spatial_solution(struct hcs *hc, struct sh_lms *sol,
+void hc_print_spatial_solution(struct hcs *hc,
+ struct sh_lms *sol,
float *sol_x, char *name,
char *dfilename,
int sol_mode, hc_boolean binary,
@@ -122,16 +131,21 @@
*/
sh_compute_spatial_basis(sol, file_dummy, FALSE,flt_dummy, &xy,
1,verbose);
- /*
- compute the scaling factors
- */
- hc_compute_solution_scaling_factors(hc,sol_mode,fac);
+
/* depth file */
dout = hc_open(dfilename,"w","hc_print_spatial_solution");
if(verbose >= 2)
fprintf(stderr,"hc_print_spatial_solution: writing depth levels to %s\n",
dfilename);
for(i=0;i < nradp2;i++){
+ /*
+
+ compute the scaling factors, those do depend on radius
+ in the case of the stresses, so leave inside loop!
+
+ */
+ hc_compute_solution_scaling_factors(hc,sol_mode,hc->r[i],fac);
+
/* write depth in [km] to dout file */
fprintf(dout,"%g\n",HC_Z_DEPTH(hc->r[i]));
for(k=0;k < 3;k++) /* pointers */
@@ -141,8 +155,10 @@
format:
- lon lat vr vt vp
+ lon lat vr vt vp OR
+ lon lat srr srt srp
+
*/
@@ -266,17 +282,21 @@
}
/*
- compute the r, theta, phi fac[3] scaling factors for the solution output
+ compute the r, theta, phi fac[3] scaling factors for the solution
+ output
*/
-void hc_compute_solution_scaling_factors(struct hcs *hc,int sol_mode,HC_PREC *fac)
+void hc_compute_solution_scaling_factors(struct hcs *hc,
+ int sol_mode,
+ HC_PREC radius,
+ HC_PREC *fac)
{
switch(sol_mode){
case HC_VEL:
fac[0]=fac[1]=fac[2] = hc->vel_scale; /* go to cm/yr */
break;
- case HC_STRESS:
- fac[0]=fac[1]=fac[2] = 1.0; /* go to ??? */
+ case HC_TRACTIONS:
+ fac[0]=fac[1]=fac[2] = hc->stress_scale/radius; /* go to MPa */
break;
default:
HC_ERROR("hc_print_spectral_solution","mode undefined");
Modified: mc/1D/hc/trunk/hc_polsol.c
===================================================================
--- mc/1D/hc/trunk/hc_polsol.c 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/hc_polsol.c 2006-07-31 21:02:14 UTC (rev 4124)
@@ -79,9 +79,7 @@
*/
hc_boolean verbose)
{
- //
- // PROGRAM POLNOV
- //
+
// ****************************************************************
// * THIS PROGRAM IS TO USED CALCULATE AND OUTPUT THE POLOIDAL *
// * COMPONENTS OF THE FLOW WITH PLATE MOTIONS *
Modified: mc/1D/hc/trunk/hc_solve.c
===================================================================
--- mc/1D/hc/trunk/hc_solve.c 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/hc_solve.c 2006-07-31 21:02:14 UTC (rev 4124)
@@ -142,11 +142,9 @@
if(verbose)
fprintf(stderr,"hc_solve: computing solution for velocities\n");
break;
- case HC_STRESS:
+ case HC_TRACTIONS:
if(verbose)
- fprintf(stderr,"hc_solve: computing solution for stresses\n");
- if(compute_geoid)
- HC_ERROR("hc_solve","cannot compute stresses and geoid");
+ fprintf(stderr,"hc_solve: computing solution for tractions\n");
break;
default:
fprintf(stderr,"hc_solve: error: solution mode %i undefined\n",
@@ -192,6 +190,11 @@
pol_sol[6*nradp2]: y1...y6 (six) poloidal solutions for each layer
tor_sol[2*nradp2]: y9 and y10 (two) toroidal solutions for each layer
+
+THESE SOLUTIONS WILL NEED TO BE SCALED WITH CONSTANTS AS GIVEN IN
+
+hc_compute_solution_scaling_factors
+
*/
void hc_sum(struct hcs *hc,
int nrad,struct sh_lms *pol_sol, struct sh_lms *tor_sol,
@@ -210,20 +213,26 @@
solution
*/
- if (solve_mode == HC_VEL){
+ switch(solve_mode){
+ case HC_VEL:
//
// velocity output requested
//
irchoose = 0; // y1 for radial
ipchoose = 1; // y2 for poloidal
itchoose = 0; // y9 for toroidal
- }else{
+ break;
+ case HC_TRACTIONS:
//
// srr srt srp stress output requested
//
irchoose = 2;// y3 for radial
ipchoose = 3;// y4 for poloidal
itchoose = 1;// y10 for toroidal
+ break;
+ default:
+ HC_ERROR("hc_sum","solve mode undefined");
+ break;
}
/*
Modified: mc/1D/hc/trunk/main.c
===================================================================
--- mc/1D/hc/trunk/main.c 2006-07-29 00:08:24 UTC (rev 4123)
+++ mc/1D/hc/trunk/main.c 2006-07-31 21:02:14 UTC (rev 4124)
@@ -17,7 +17,34 @@
$Id: main.c,v 1.13 2006/01/22 01:11:34 becker Exp becker $
+>>> SOME COMMENTS FROM THE ORIGINAL CODE <<<
+C * It uses the following Numerical Recipes (Press et al.) routines:
+C four1, realft, gauleg, rk4, rkdumb, ludcmp, lubksb;
+C !!!!!!!!!!!!!!!!!!! rkqc, odeint !!!!!!!!!! take out !!!!!!
+C and the following routines by R.J. O'Connell:
+C shc2d, shd2c, ab2cs, cs2ab, plmbar, plvel2sh, pltgrid, pltvel,
+C vshd2c, plmbar1.
+C Further subroutines are: kinsub, evalpa,
+C torsol (all based on "kinflow" by Hager & O'Connell),
+C densub and evppot (based on "denflow" by Hager & O'Connell),
+C sumsub (based on "sumkd" by Hager & O'Connell, but
+C substantially speeded up),
+C convert, derivs and shc2dd (based on R.J. O'Connell's shc2d).
+C
+C bugs found:
+C * The combination of (1) high viscosity lithosphere
+C (2) compressible flow (3) kinematic (plate driven) flow
+C doesn't work properly. The problem presumably only occurs
+C at degree 1 (I didn't make this sure) but this is sufficient
+C to screw up everything. It will usually work to reduce the
+C contrast between lithospheric and asthenospheric viscosity.
+C Then make sure that (1) two viscosity structures give similar
+C results for incompressible and (2) incompressible and compressible
+C reduced viscosity look similar (e.g. anomalous mass flux vs. depth)
+
+<<< END OLD COMMENTS
+
*/
int main(int argc, char **argv)
@@ -27,9 +54,8 @@
struct sh_lms *sol_spectral=NULL, *geoid = NULL; /* solution expansions */
int nsol,lmax;
FILE *out;
- hc_boolean compute_geoid = TRUE; /* compute the geoid (only works for velocity solution) */
struct hc_parameters p[1]; /* parameters */
- char filename[HC_CHAR_LENGTH];
+ char filename[HC_CHAR_LENGTH],file_prefix[10];
float *sol_spatial = NULL; /* spatial solution,
e.g. velocities */
/*
@@ -83,9 +109,11 @@
*/
hc_init_main(model,SH_RICK,p);
nsol = (model->nrad+2) * 3; /* number of solution (r,pol,tor)*(nlayer+2) */
- if(p->free_slip)
+ if(p->free_slip) /* maximum degree is determined by the
+ density expansion */
lmax = model->dens_anom[0].lmax;
- else
+ else /* max degree is determined by the
+ plate velocities */
lmax = model->pvel[0].lmax; /* shouldn't be larger than that*/
@@ -103,7 +131,7 @@
*/
sh_allocate_and_init(&sol_spectral,nsol,lmax,model->sh_type,1,
p->verbose);
- if(compute_geoid)
+ if(p->compute_geoid)
/* make room for geoid solution */
sh_allocate_and_init(&geoid,1,model->dens_anom[0].lmax,
model->sh_type,0,p->verbose);
@@ -112,7 +140,7 @@
solve poloidal and toroidal part and sum
*/
hc_solve(model,p->free_slip,p->solution_mode,sol_spectral,
- TRUE,TRUE,TRUE,p->print_pt_sol,compute_geoid,geoid,
+ TRUE,TRUE,TRUE,p->print_pt_sol,p->compute_geoid,geoid,
p->verbose);
/*
@@ -122,10 +150,18 @@
/*
output of spherical harmonics solution
*/
+ switch(p->solution_mode){
+ case HC_VEL:
+ sprintf(file_prefix,"vel");break;
+ case HC_TRACTIONS:
+ sprintf(file_prefix,"str");break;
+ default:
+ HC_ERROR(argv[0],"solution mode undefined");break;
+ }
if(p->sol_binary_out)
- sprintf(filename,"%s",HC_SOLOUT_FILE_BINARY);
+ sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_BINARY);
else
- sprintf(filename,"%s",HC_SOLOUT_FILE_ASCII);
+ sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_ASCII);
if(p->verbose)
fprintf(stderr,"%s: writing spherical harmonics solution to %s\n",
argv[0],filename);
@@ -134,7 +170,7 @@
p->solution_mode,
p->sol_binary_out,p->verbose);
fclose(out);
- if(compute_geoid){
+ if(p->compute_geoid){
/*
print geoid solution
*/
@@ -159,10 +195,10 @@
/*
output of spatial solution
*/
+ sprintf(filename,"%s.%s",file_prefix,HC_SPATIAL_SOLOUT_FILE);
/* print lon lat z v_r v_theta v_phi */
hc_print_spatial_solution(model,sol_spectral,sol_spatial,
- HC_SPATIAL_SOLOUT_FILE,
- HC_LAYER_OUT_FILE,
+ filename,HC_LAYER_OUT_FILE,
p->solution_mode,p->sol_binary_out,
p->verbose);
}
@@ -172,7 +208,7 @@
*/
sh_free_expansion(sol_spectral,nsol);
- if(compute_geoid)
+ if(p->compute_geoid)
sh_free_expansion(geoid,1);
free(sol_spatial);
More information about the cig-commits
mailing list