[cig-commits] r12213 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Fri Jun 13 17:44:32 PDT 2008
Author: becker
Date: 2008-06-13 17:44:32 -0700 (Fri, 13 Jun 2008)
New Revision: 12213
Added:
mc/1D/hc/trunk/hc_constants.h
Modified:
mc/1D/hc/trunk/Makefile
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/main.c
Log:
Moved several more constants into a single file to make it easier to work out
the Earth model parameters.
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2008-06-14 00:21:26 UTC (rev 12212)
+++ mc/1D/hc/trunk/Makefile 2008-06-14 00:44:32 UTC (rev 12213)
@@ -25,7 +25,7 @@
BDIR = ../hc/bin/$(ARCH)/
# include files
-OINCS = hc.h hc_filenames.h sh.h
+OINCS = hc.h hc_filenames.h sh.h hc_constants.h
#
# Healpix stuff, comment out all if not wanted
#
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2008-06-14 00:21:26 UTC (rev 12212)
+++ mc/1D/hc/trunk/hc.h 2008-06-14 00:44:32 UTC (rev 12213)
@@ -68,8 +68,14 @@
*/
#include "prem.h"
+/*
+Other Earth model, and modeling constants
+*/
+#include "hc_constants.h"
+
+
/*
dealing with velocity grids
@@ -429,16 +435,8 @@
#endif
/*
-other constants
*/
-// now taken from earth model
-#define HC_RE_KM 6371.0 /* radius(Earth) in [km] */
-
-//#define HC_RCMB_ND 0.546225 /* non-dim radius, ~10km above
-// CMB */
-
-#define HC_TIMESCALE_YR 1e6 /* timescale [yr] */
/*
convert depth (>0) in km to non-dimensionalizd radius
*/
@@ -448,14 +446,6 @@
*/
#define HC_Z_DEPTH(x) ((HC_RE_KM * (1.0-(x))))
-/*
-default constant scaling for input density file, use 0.01 for %
-
-*/
-
-#define HC_DENSITY_SCALING 0.01
-
-
#define __HC_READ_HEADER_FILE__
#endif
Added: mc/1D/hc/trunk/hc_constants.h
===================================================================
--- mc/1D/hc/trunk/hc_constants.h (rev 0)
+++ mc/1D/hc/trunk/hc_constants.h 2008-06-14 00:44:32 UTC (rev 12213)
@@ -0,0 +1,45 @@
+/*
+
+Earth model constants
+
+
+*/
+// now taken from earth model
+#define HC_RE_KM 6371.0 /* radius(Earth) in [km] */
+
+
+#define HC_TIMESCALE_YR 1e6 /* timescale [yr] */
+
+
+#define HC_VISNOR 1e21 /* reference viscosity [Pas] */
+#define HC_GACC 10.0e2 /* gravitational acceleration [cm/s2] */
+#define HC_CAPITAL_G 6.6742e-11 /* gravitational constant [Nm2/kg2]*/
+
+#define HC_SECYR 3.1556926e7 /* seconds/year */
+
+#define HC_AVG_DEN_MANTLE 4.4488 /* in g/cm^3 */
+#define HC_AVG_DEN_CORE 11.60101
+
+
+/*
+
+other modeling constants or initial values
+
+*/
+/*
+
+default constant scaling for input density file, use 0.01 for %
+
+*/
+
+#define HC_DENSITY_SCALING 0.01
+
+/* default dln v_s/dln rho density scaling */
+
+#define HC_D_LOG_V_D_LOG_D 0.2
+
+
+/* length scale factor for velocity I/O, in units of m
+ use 0.01, if input/output is to be in cm/yr
+*/
+#define HC_VEL_IO_SCALE 0.01
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2008-06-14 00:21:26 UTC (rev 12212)
+++ mc/1D/hc/trunk/hc_init.c 2008-06-14 00:44:32 UTC (rev 12213)
@@ -28,9 +28,9 @@
p->compute_geoid_correlations = FALSE; /* compute the geoid
correlation with
refernece only */
- p->dens_anom_scale = 0.2; /* default density anomaly scaling to
- go from PREM percent traveltime
- anomalies to density anomalies */
+ p->dens_anom_scale = HC_D_LOG_V_D_LOG_D ; /* default density anomaly scaling to
+ go from PREM percent traveltime
+ anomalies to density anomalies */
p->read_short_dens_sh = FALSE; /* read the density anomaly file in
the short format of Becker &
Boschi (2002)? */
@@ -252,9 +252,9 @@
constants
*/
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]*/
+ hc->visnor = HC_VISNOR; /* normalizing viscosity [Pas]*/
+ hc->gacc = HC_GACC; /* gravitational acceleration [cm/s2]*/
+ hc->g = HC_CAPITAL_G; /* gravitational constant [Nm2/kg2]*/
/*
@@ -264,12 +264,14 @@
hc->re = hc->prem->r0;
if(fabs(hc->re - (HC_RE_KM * 1e3)) > 1e-7)
HC_ERROR("hc_init_constants","Earth radius mismatch")
- hc->secyr = 3.1556926e7; /* seconds/year */
+
+ hc->secyr = HC_SECYR; /* seconds/year */
+
/*
those are in g/cm^3
*/
- hc->avg_den_mantle = 4.4488;
- hc->avg_den_core = 11.60101;
+ hc->avg_den_mantle = HC_AVG_DEN_MANTLE;
+ hc->avg_den_core = HC_AVG_DEN_CORE;
/*
take the CMB radius from the Earth model
@@ -288,7 +290,7 @@
velocity scale if input is in [cm/yr], works out to be ~0.11
*/
- hc->vel_scale = hc->re*PIOVERONEEIGHTY/hc->timesc*100;
+ hc->vel_scale = hc->re*PIOVERONEEIGHTY/hc->timesc/HC_VEL_IO_SCALE;
/*
stress scaling, will later be divided by non-dim radius, to go
@@ -296,7 +298,7 @@
*/
hc->stress_scale = (PIOVERONEEIGHTY * hc->visnor / hc->secyr)/
- (hc->timesc * 1e6);
+ (hc->timesc * HC_TIMESCALE_YR);
hc->const_init = TRUE;
@@ -468,10 +470,18 @@
*/
hc_vecrealloc(&hc->rvisc,4,"hc_assign_viscosity");
hc_vecrealloc(&hc->visc,4,"hc_assign_viscosity");
+ /* number of layers */
hc->nvis = 4;
- hc->rvisc[0] = hc->r_cmb;hc->rvisc[1] = p->rlayer[0];hc->rvisc[2] = p->rlayer[1];hc->rvisc[3] = p->rlayer[2];
- for(i=0;i < hc->nvis;i++)
+ /* radii */
+ hc->rvisc[0] = hc->r_cmb;
+ hc->rvisc[1] = p->rlayer[0];
+ hc->rvisc[2] = p->rlayer[1];
+ hc->rvisc[3] = p->rlayer[2];
+
+ for(i=0;i < hc->nvis;i++){
hc->visc[i] = elayer[i];
+ //fprintf(stderr,"%11g %11g\n",hc->rvisc[i],hc->visc[i]);
+ }
if(p->verbose)
fprintf(stderr,"hc_assign_viscosity: assigned four layer viscosity: %.2e %.2e %.2e %.2e\n",
hc->visc[0],hc->visc[1],hc->visc[2],hc->visc[3]);
Modified: mc/1D/hc/trunk/main.c
===================================================================
--- mc/1D/hc/trunk/main.c 2008-06-14 00:21:26 UTC (rev 12212)
+++ mc/1D/hc/trunk/main.c 2008-06-14 00:44:32 UTC (rev 12213)
@@ -59,7 +59,7 @@
float *sol_spatial = NULL; /* spatial solution,
e.g. velocities */
HC_PREC corr[2]; /* correlations */
- HC_PREC vl[3][3],v1,v2,v3; /* for viscosity scans */
+ HC_PREC vl[4][3],v[4],dv; /* for viscosity scans */
/*
@@ -218,35 +218,53 @@
case HC_SOLVER_MODE_VISC_SCAN: /* scan through viscosities for a
four layer problem */
- /* log bounds */
- vl[0][0]= -3;vl[0][1]=5+1e-5;vl[0][2]=.1; /* 0..100 layer log bounds and spacing */
- vl[1][0]= -3;vl[1][1]=5+1e-5;vl[1][2]=.1; /* 100..410 */
- vl[2][0]= -3;vl[2][1]=5+1e-5;vl[2][2]=.1; /* 660 ... 2871 */
- /* linear, in units of reference */
- p->elayer[1] = 1.0; /* 410..660 reference, remains unchanged */
+
+ /* parameter space log bounds */
+
+ dv = .1; /* spacing */
+
+ vl[0][0]= -3;vl[0][1]=3+1e-5;vl[0][2]=dv; /* 0..100 layer log bounds and spacing */
+ vl[1][0]= -3;vl[1][1]=3+1e-5;vl[1][2]=dv; /* 100..410 */
+ if(p->free_slip){
+ vl[2][0]= 0;vl[2][1]=0+1e-5;vl[2][2]=dv; /* for free slip,
+ only relative
+ viscosisites
+ matter for
+ correlation */
+ }else{
+ vl[2][0]= -3;vl[2][1]=3+1e-5;vl[2][2]=dv; /* need to actually
+ loop 410 .660 */
+ }
+ vl[3][0]= -3;vl[3][1]=3+1e-5;vl[3][2]=dv; /* 660 ... 2871 */
+
+
/* */
solved=0;
- for(v1=vl[0][0];v1 <= vl[0][1];v1 += vl[0][2])
- for(v2=vl[1][0];v2 <= vl[1][1];v2 += vl[1][2])
- for(v3=vl[2][0];v3 <= vl[2][1];v3 += vl[2][2]){
- /* layer viscosity structure */
- p->elayer[0] = pow(10,v3);p->elayer[2] = pow(10,v2);p->elayer[3]=pow(10,v1);
- hc_assign_viscosity(model,HC_INIT_E_FOUR_LAYERS,p->elayer,p);
- /* print viscosities of 0...100, 100...410, and 660...2871 layer */
- fprintf(stdout,"%14.7e %14.7e %14.7e\t",p->elayer[3],p->elayer[2],p->elayer[0]);
- /* compute solution */
- hc_solve(model,p->free_slip,p->solution_mode,sol_spectral,
- (solved)?(FALSE):(TRUE), /* density changed? */
- (solved)?(FALSE):(TRUE), /* plate velocity changed? */
- TRUE, /* viscosity changed */
- FALSE,p->compute_geoid,geoid,
- p->verbose);
- /* only output are the geoid correlations, for now */
- hc_compute_correlation(geoid,p->ref_geoid,corr,1,p->verbose);
- fprintf(stdout,"%10.7f %10.7f ",corr[0],corr[1]);
- fprintf(stdout,"\n");
- solved++;
- }
+ for(v[0]=vl[0][0];v[0] <= vl[0][1];v[0] += vl[0][2])
+ for(v[1]=vl[1][0];v[1] <= vl[1][1];v[1] += vl[1][2])
+ for(v[2]=vl[2][0];v[2] <= vl[2][1];v[2] += vl[2][2])
+ for(v[3]=vl[3][0];v[3] <= vl[3][1];v[3] += vl[3][2]){
+ /* layer viscosity structure */
+ p->elayer[0] = pow(10,v[3]); /* bottom */
+ p->elayer[1] = pow(10,v[2]); /* 660..410 */
+ p->elayer[2] = pow(10,v[1]); /* 410..100 */
+ p->elayer[3] = pow(10,v[0]); /* 100..0 */
+ hc_assign_viscosity(model,HC_INIT_E_FOUR_LAYERS,p->elayer,p);
+ /* print viscosities of 0...100, 100...410, 410 ... 660 and 660...2871 layer */
+ fprintf(stdout,"%14.7e %14.7e %14.7e %14.7e\t",p->elayer[3],p->elayer[2],p->elayer[1],p->elayer[0]);
+ /* compute solution */
+ hc_solve(model,p->free_slip,p->solution_mode,sol_spectral,
+ (solved)?(FALSE):(TRUE), /* density changed? */
+ (solved)?(FALSE):(TRUE), /* plate velocity changed? */
+ TRUE, /* viscosity changed */
+ FALSE,p->compute_geoid,geoid,
+ p->verbose);
+ /* only output are the geoid correlations, for now */
+ hc_compute_correlation(geoid,p->ref_geoid,corr,1,p->verbose);
+ fprintf(stdout,"%10.7f %10.7f ",corr[0],corr[1]);
+ fprintf(stdout,"\n");
+ solved++;
+ }
break;
default:
HC_ERROR("hc","solver mode undefined");
More information about the cig-commits
mailing list