[cig-commits] r11868 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Tue Apr 29 12:34:22 PDT 2008
Author: becker
Date: 2008-04-29 12:34:22 -0700 (Tue, 29 Apr 2008)
New Revision: 11868
Added:
mc/1D/hc/trunk/compute_and_extract_tractions
Modified:
mc/1D/hc/trunk/ggrd_grdtrack_util.c
mc/1D/hc/trunk/ggrd_struc.h
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_polsol.c
mc/1D/hc/trunk/rotvec2vel.c
mc/1D/hc/trunk/sh_syn.c
Log:
Fixed minor typos, added another example script
Added: mc/1D/hc/trunk/compute_and_extract_tractions
===================================================================
--- mc/1D/hc/trunk/compute_and_extract_tractions (rev 0)
+++ mc/1D/hc/trunk/compute_and_extract_tractions 2008-04-29 19:34:22 UTC (rev 11868)
@@ -0,0 +1,52 @@
+#!/bin/bash
+
+#
+# simple example on how to compute radial tractions and extract them from the hc solution
+#
+# data dir
+ddir=$HOME/progs/src/seatree/py-drivers/py-hc/example_data/
+#
+# viscosity structure
+vfile=$ddir/viscosity/visc.C
+#
+# density model
+dfile=$ddir/tomography/smean.31.m.ab
+dformat="-dshs" # short SH format as in Becker & Boschi
+dscale=0.25 # dln rho/dln v_S factor
+#
+# plate velocities
+pfile=$ddir/pvelocity/nnr_nuvel1a.smoothed.sh.dat
+
+# verbosity
+verbose="-vvv"
+
+# options: tractions, and no geoid outut
+opt="-rtrac -ng"
+
+# free slip
+hc $verbose -dens $dfile $dformat -ds $dscale -vf $vfile $opt
+mv rtrac.sol.bin trac.fs.sol
+
+# plates
+hc $verbose -dens $dfile $dformat -ds $dscale -vf $vfile -pvel $pfile $opt
+mv rtrac.sol.bin trac.p.sol
+
+# extract
+layer=36
+for t in p fs;do
+ hc_extract_sh_layer trac.$t.sol -2 4 > tmp.$$
+ z=`gawk '{if(NR==l)print($2)}' l=$layer tmp.$$`
+ rm tmp.$$
+ echo solution trac.$t.sol layer $layer depth $z
+ # radial
+ hc_extract_sh_layer trac.$t.sol $layer 1 | sh_syn 0 0 > tmp.$$.r
+ # theta phi
+ hc_extract_sh_layer trac.$t.sol $layer 2 | sh_syn 0 1 > tmp.$$.tp
+
+ # combine to lon lat srr srt srp format
+ echo $0: writing tractions to trac.$t.$z.dat
+ paste tmp.$$.r tmp.$$.tp | gawk '{print($1,$2,$3,$6,$7)}' > trac.$t.$z.dat
+
+ rm tmp.$$.r tmp.$$.tp
+done
+
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2008-04-29 18:02:53 UTC (rev 11867)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2008-04-29 19:34:22 UTC (rev 11868)
@@ -24,6 +24,7 @@
void ggrd_init_master(struct ggrd_master *ggrd)
{
ggrd->mat_control = ggrd->mat_control_init = 0;
+ ggrd->ray_control = ggrd->ray_control_init = 0;
ggrd->vtop_control = ggrd->vtop_control_init = 0;
ggrd->age_control = ggrd->age_control_init = 0;
ggrd->nage = 0;ggrd->age_bandlim = 200.;
Modified: mc/1D/hc/trunk/ggrd_struc.h
===================================================================
--- mc/1D/hc/trunk/ggrd_struc.h 2008-04-29 18:02:53 UTC (rev 11867)
+++ mc/1D/hc/trunk/ggrd_struc.h 2008-04-29 19:34:22 UTC (rev 11868)
@@ -131,15 +131,20 @@
struct ggrd_master{ /* master structure */
int mat_control,mat_control_init;
+ int ray_control,ray_control_init;
int vtop_control,vtop_control_init;
int age_control,age_control_init;
char mat_file[1000];
+ char ray_file[1000];
char vtop_dir[1000];
char age_file[1000];
/* grid structures */
struct ggrd_gt *mat; /* material grids */
+ /* rayleigh number grids */
+ struct ggrd_gt *ray;
+
/* surface velocity */
struct ggrd_gt *svp,*svt; /* phi/theta surface velocities */
/* age stuff */
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2008-04-29 18:02:53 UTC (rev 11867)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2008-04-29 19:34:22 UTC (rev 11868)
@@ -130,7 +130,7 @@
void rick_plmbar1(float *, float *, int, int, float, struct rick_module *);
void rick_gauleg(float, float, float *, float *, int);
/* rotvec2vel.c */
-FILE *myopen(const char *, const char *);
+FILE *rv_myopen(const char *, const char *);
/* sh_ana.c */
/* shana_sh.c */
void shana_compute_allplm(int, int, double *, double *, struct shana_module *);
Modified: mc/1D/hc/trunk/hc_polsol.c
===================================================================
--- mc/1D/hc/trunk/hc_polsol.c 2008-04-29 18:02:53 UTC (rev 11867)
+++ mc/1D/hc/trunk/hc_polsol.c 2008-04-29 19:34:22 UTC (rev 11868)
@@ -340,6 +340,7 @@
// APPEND A FINAL RVISC_LOCAL = 1.0 TO PREVENT OUT OF BOUNDS
//
hc->rvisc[hc->nvis] = 1.0;
+
//
// INITIALIZE INDEX,IVIS,NIH
//
@@ -360,7 +361,9 @@
<HC_EPS_PREC)){
hit = TRUE; /* bailout here */
}else{
- /* normal operation */
+ /*
+ normal operation
+ */
hc->pvisc[hc->nprops] = hc->visc[ivis];
hc->den[hc->nprops] = 0.0;
hc->qwrite[hc->nprops] = FALSE;
@@ -415,6 +418,8 @@
}
} /* end of nrad loop */
hc->den[hc->nprops] = 0.0;
+ hc->pvisc[hc->nprops] = hc->pvisc[hc->nprops-1]; /* to look nicer */
+
/*
number of propagators is now nprops+1
@@ -486,6 +491,8 @@
*/
}
hc->rprops[hc->nprops+1] = 1.0;
+
+
if(verbose >= 3)
for(i=0;i < hc->nprops+2;i++)
fprintf(stderr,"i: %3i nprops: %3i r(i): %11g rho: %11g\n",
Modified: mc/1D/hc/trunk/rotvec2vel.c
===================================================================
--- mc/1D/hc/trunk/rotvec2vel.c 2008-04-29 18:02:53 UTC (rev 11867)
+++ mc/1D/hc/trunk/rotvec2vel.c 2008-04-29 19:34:22 UTC (rev 11868)
@@ -28,7 +28,7 @@
//#define VELFACTOR (6371009.0/1.0e6*1.0e2)
//#define PIOVERONEEIGHTY_TIMES_VELFACTOR (PIOVERONEEIGHTY*VELFACTOR)
#define PIOVERONEEIGHTY_TIMES_VELFACTOR 11.1195083724191
-FILE *myopen(const char *,const char *);
+FILE *rv_myopen(const char *,const char *);
int main(int argc, char *argv[])
@@ -71,7 +71,7 @@
for(i=0;i<CODED_PLATES;i++)
assigned[i]=0;
- vel=myopen(argv[1],"r");
+ vel=rv_myopen(argv[1],"r");
rotvel=(double *)malloc(sizeof(double)*3*CODED_PLATES);
stats=(int *)calloc(CODED_PLATES,sizeof(int));
@@ -275,7 +275,7 @@
fprintf(stderr,"%s: done\n",argv[0]);
return 0;
}
-FILE *myopen(const char *name,const char *modus)
+FILE *rv_myopen(const char *name,const char *modus)
{
FILE *tmp;
if( (tmp = (FILE *)fopen(name,modus)) == NULL)
Modified: mc/1D/hc/trunk/sh_syn.c
===================================================================
--- mc/1D/hc/trunk/sh_syn.c 2008-04-29 18:02:53 UTC (rev 11867)
+++ mc/1D/hc/trunk/sh_syn.c 2008-04-29 19:34:22 UTC (rev 11868)
@@ -63,7 +63,7 @@
sscanf(argv[8],"%lf",&dy);
else
dy = dx;
- if((argc > 8)|| (argc < 0)){
+ if((argc > 9)|| (argc < 0)){
fprintf(stderr,"usage: %s [short_format, %i] [short_ivec, %i] [w, %g] [e, %g] [s, %g] [n, %g] [dx, %g] [dy, dx] (in that order)\n",
argv[0],short_format,short_format_ivec,w,e,s,n,dx);
fprintf(stderr,"short_format:\n\t0: expects regular format with long header\n");
More information about the cig-commits
mailing list