[cig-commits] r18781 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Sun Jul 17 23:12:08 PDT 2011
Author: becker
Date: 2011-07-17 23:12:07 -0700 (Sun, 17 Jul 2011)
New Revision: 18781
Modified:
mc/1D/hc/trunk/Makefile
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_extract_sh_layer.c
mc/1D/hc/trunk/hc_extract_spatial.c
mc/1D/hc/trunk/hc_filenames.h
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/hc_input.c
mc/1D/hc/trunk/hc_misc.c
mc/1D/hc/trunk/hc_output.c
mc/1D/hc/trunk/main.c
mc/1D/hc/trunk/sh.h
mc/1D/hc/trunk/sh_exp.c
Log:
implemented proper vtk output and fixed a bug
in copy routine for spherical harmonics, apparently
never caused a problem.
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2011-07-17 13:59:04 UTC (rev 18780)
+++ mc/1D/hc/trunk/Makefile 2011-07-18 06:12:07 UTC (rev 18781)
@@ -1,3 +1,4 @@
+#CFLAGS = -g
#
#
# makefile for experimental Hager & O'Connell routines
@@ -11,7 +12,6 @@
#
# EDIT HERE FOR GMT VERSION
#
-
include Makefile.include
#
#
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2011-07-17 13:59:04 UTC (rev 18780)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2011-07-18 06:12:07 UTC (rev 18781)
@@ -55,13 +55,14 @@
void hc_read_geoid(struct hc_parameters *);
void hc_select_pvel(double, struct pvels *, struct sh_lms *, unsigned short);
/* hc_input.c */
-void hc_read_sh_solution(struct hcs *, struct sh_lms **, FILE *, unsigned short, unsigned short);
+int hc_read_sh_solution(struct hcs *, struct sh_lms **, FILE *, unsigned short, unsigned short);
/* hc_matrix.c */
void hc_ludcmp_3x3(double [3][3], int, int *);
void hc_lubksb_3x3(double [3][3], int, int *, double *);
/* hc_misc.c */
void hc_dvecalloc(double **, int, char *);
void hc_svecalloc(float **, int, char *);
+void hc_ivecalloc(int **, int, char *);
void hc_vecalloc(double **, int, char *);
void hc_scmplx_vecalloc(struct hc_scmplx **, int, char *);
void hc_svecrealloc(float **, int, char *);
@@ -83,6 +84,7 @@
void lonlatpv2cv(float, float, float *, float *);
void lonlatpv2cv_with_base(float *, double *, float *);
void calc_polar_base_at_theta_phi(float, float, double *);
+void hc_linear_interpolate(double *, int, double, int *, int *, double *, double *);
/* hc_output.c */
void hc_print_spectral_solution(struct hcs *, struct sh_lms *, FILE *, int, unsigned short, unsigned short);
void hc_print_sh_scalar_field(struct sh_lms *, FILE *, unsigned short, unsigned short, unsigned short);
@@ -97,11 +99,13 @@
void hc_compute_solution_scaling_factors(struct hcs *, int, double, 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);
-void hc_print_vtk(FILE *, float *, float *, int, int, unsigned short);
+void hc_print_vtk(FILE *, float *, float *, int, int, unsigned short, int, float *, int, int);
void hc_print_be_float(float *, int, FILE *, unsigned short);
+void hc_print_be_int(int *, int, FILE *, unsigned short);
unsigned short hc_is_little_endian(void);
void hc_flip_byte_order(void *, size_t);
void hc_flipit(void *, void *, size_t);
+void hc_print_dens_anom(struct hcs *, FILE *, unsigned short, unsigned short);
/* hc_polsol.c */
void hc_polsol(struct hcs *, int, double *, int, double *, unsigned short, struct sh_lms *, unsigned short, int, double *, double *, unsigned short, struct sh_lms *, struct sh_lms *, unsigned short, struct sh_lms *, unsigned short, unsigned short);
/* hc_propagator.c */
@@ -196,6 +200,7 @@
void sh_add_coeff(struct sh_lms *, int, int, int, unsigned short, double *);
void sh_copy_lms(struct sh_lms *, struct sh_lms *);
void sh_aexp_equals_bexp_coeff(struct sh_lms *, struct sh_lms *);
+void sh_c_is_a_plus_b_coeff(struct sh_lms *, struct sh_lms *, struct sh_lms *);
void sh_scale_expansion_l_factor(struct sh_lms *, double *);
void sh_scale_expansion(struct sh_lms *, double);
/* sh_extract_layer.c */
Modified: mc/1D/hc/trunk/hc_extract_sh_layer.c
===================================================================
--- mc/1D/hc/trunk/hc_extract_sh_layer.c 2011-07-17 13:59:04 UTC (rev 18780)
+++ mc/1D/hc/trunk/hc_extract_sh_layer.c 2011-07-18 06:12:07 UTC (rev 18781)
@@ -12,12 +12,12 @@
int main(int argc, char **argv)
{
- int ilayer,nsol,i,mode,shps,nset=1,loop,i1,i2;
+ int ilayer,nsol,i,mode,shps=1,nset=1,loop,i1,i2,shps_read;
FILE *in;
struct sh_lms *sol=NULL;
struct hcs *model;
HC_PREC fac[3] = {1.0,1.0,1.0};
- hc_boolean binary = TRUE, verbose = FALSE, short_format = FALSE;
+ hc_boolean binary = TRUE, verbose = TRUE, short_format = FALSE;
hc_struc_init(&model);
/*
deal with parameters
@@ -65,9 +65,9 @@
read in solution
*/
in = ggrd_open(argv[1],"r","hc_extract_sh_layer");
- hc_read_sh_solution(model,&sol,in,binary,verbose);
+ shps_read = hc_read_sh_solution(model,&sol,in,binary,verbose);
fclose(in);
- nsol = model->nradp2 * 3;
+ nsol = model->nradp2 * shps_read;
/*
deal with selection
*/
@@ -90,13 +90,22 @@
}else{
i1=ilayer-1;i2 = i1;
}
- /* detect number of expansions */
+ /*
+ detect number of expansions
+ */
if((mode == 1)||(mode == 5)||(mode == 6))
shps = 1;
else if(mode == 2)
shps = 2;
else if(mode == 3)
shps = 3;
+
+ if(shps > shps_read){
+ fprintf(stderr,"%s: solution file only had %i expansions, mode %i requests %i\n",
+ argv[0],shps_read,mode,shps);
+ exit(-1);
+
+ }
for(ilayer=i1;ilayer <= i2;ilayer++){
/*
output
@@ -105,7 +114,7 @@
/* SH header */
if(short_format && loop)
fprintf(stdout,"%g\n",HC_Z_DEPTH(model->r[ilayer]));
- sh_print_parameters_to_file((sol+ilayer*3),shps,
+ sh_print_parameters_to_file((sol+ilayer*shps_read),shps,
ilayer,nset,(HC_PREC)(HC_Z_DEPTH(model->r[ilayer])),
stdout,short_format,FALSE,verbose);
}
@@ -115,21 +124,21 @@
if(verbose)
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);
+ sh_print_coefficients_to_file((sol+ilayer*shps_read),shps,stdout,fac,FALSE,verbose);
break;
case 2:
/* */
if(verbose)
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);
+ sh_print_coefficients_to_file((sol+ilayer*shps_read+1),shps,stdout,fac,FALSE,verbose);
break;
case 3:
/* mode == 3 */
if(verbose)
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);
+ sh_print_coefficients_to_file((sol+ilayer*shps_read),shps,stdout,fac,FALSE,verbose);
break;
case 4:
fprintf(stdout,"%5i %11g\n",ilayer,HC_Z_DEPTH(model->r[ilayer]));
@@ -139,14 +148,14 @@
if(verbose)
fprintf(stderr,"%s: printing x_pol 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);
+ sh_print_coefficients_to_file((sol+ilayer*shps_read+1),shps,stdout,fac,FALSE,verbose);
break;
case 6:
/* */
if(verbose)
fprintf(stderr,"%s: printing 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+2),shps,stdout,fac,FALSE,verbose);
+ sh_print_coefficients_to_file((sol+ilayer*shps_read+2),shps,stdout,fac,FALSE,verbose);
break;
default:
Modified: mc/1D/hc/trunk/hc_extract_spatial.c
===================================================================
--- mc/1D/hc/trunk/hc_extract_spatial.c 2011-07-17 13:59:04 UTC (rev 18780)
+++ mc/1D/hc/trunk/hc_extract_spatial.c 2011-07-18 06:12:07 UTC (rev 18781)
@@ -5,20 +5,20 @@
extract part of a solution of a HC run and convert to spatial
-
*/
int main(int argc, char **argv)
{
- int ilayer,nsol,mode,shps,loop,i1,i2,ivec,lc,ndata,npoints,i,j,
- poff,ilayer3;
+ int ilayer,nvsol,ndsol=0,mode,shps,loop,i1,i2,nlat,nlon,
+ ivec,lc,ndata,ndata_all,ndata_d,npoints,i,j,
+ poff,shps_read=0,shps_read_d=0;
FILE *in;
- struct sh_lms *sol=NULL;
+ struct sh_lms *vsol=NULL,*dsol=NULL;
struct hcs *model;
HC_PREC zlabel;
- hc_boolean binary_in = TRUE, verbose = FALSE;
- float *data,*plm=NULL,*xpos,*xvec,lon,lat,theta,phi,xtmp[3],
- pvec[3];
+ hc_boolean binary_in = TRUE, verbose = FALSE,read_dsol=FALSE;
+ float *data,*plm=NULL,*xpos,*xvec,lon,lat,theta,phi,xtmp[3],pvec[3],
+ *xscalar;
double polar_base[9];
hc_struc_init(&model);
/*
@@ -34,10 +34,16 @@
sscanf(argv[2],"%i",&ilayer);
sscanf(argv[3],"%i",&mode);
break;
+ case 5:
+ sscanf(argv[2],"%i",&ilayer);
+ sscanf(argv[3],"%i",&mode);
+ read_dsol = TRUE;
+ break;
default:
- fprintf(stderr,"%s: usage\n%s sol.file layer [mode,%i] \n\n",
+ fprintf(stderr,"%s: usage\n%s sol.file layer [mode,%i] [scalar.sol]\n\n",
argv[0],argv[0],mode);
- fprintf(stderr,"extracts spatial solution x (velocity or stress) from HC run\n");
+ fprintf(stderr,"extracts spatial solution (velocity or stress, v) from output file sol.file\n");
+ fprintf(stderr," if scalar.sol argument is given, will also read in a scalar for VTK output\n");
fprintf(stderr,"layer: 1...nset\n");
fprintf(stderr,"\tif ilayer= 1..nset, will print one layer\n");
fprintf(stderr,"\t -1, will select nset (the top layer)\n");
@@ -47,20 +53,20 @@
fprintf(stderr,"\t 2, will print lon lat z v_theta v_phi \n");
fprintf(stderr,"\t 3, will print lon lat z v_r v_theta v_phi\n");
fprintf(stderr,"\t 4, will print the depth levels of all layers\n");
- fprintf(stderr,"\t 5, compute all depth levels and write VTK file, ASCII\n");
- fprintf(stderr,"\t 6, compute all depth levels and write VTK file, BINARY\n");
+ fprintf(stderr,"\t 5, compute all depth levels (set ilayer=-2) and write VTK file, ASCII\n");
+ fprintf(stderr,"\t 6, compute all depth levels (set ilayer=-2) and write VTK file, BINARY\n");
exit(-1);
break;
}
if((mode == 4)||(mode==5)||(mode==6))
ilayer = -2;
/*
- read in solution
+ read in velocity/traction solution
*/
- in = ggrd_open(argv[1],"r","hc_extract_sh_layer");
- hc_read_sh_solution(model,&sol,in,binary_in,verbose);
+ in = ggrd_open(argv[1],"r","hc_extract_spatial");
+ shps_read = hc_read_sh_solution(model,&vsol,in,binary_in,verbose);
fclose(in);
- nsol = model->nradp2 * 3;
+ nvsol = model->nradp2 * shps_read;
/*
deal with selection
*/
@@ -92,14 +98,44 @@
}else{
shps = 1;
}
- /* room for spatial expansion */
- npoints = (sol+i1*3)->npoints;
- ndata = npoints * shps;
- if((mode == 5)||(mode==6)) /* save all layers */
- hc_svecalloc(&data,model->nradp2 * ndata,"hc_extract_spatial");
- else
- hc_svecalloc(&data,ndata,"hc_extract_spatial");
- for(lc=0,ilayer=i1,ilayer3=ilayer*3;ilayer <= i2;ilayer++,lc++,ilayer3+=3){
+ if(shps > shps_read){
+ fprintf(stderr,"%s: solution file only had %i expansions, mode %i requests %i\n",
+ argv[0],shps_read,mode,shps);
+ exit(-1);
+ }
+ /*
+ density solution or other scalar
+ */
+ if(read_dsol){
+ if((mode != 5)&&(mode != 6))
+ HC_ERROR("hc_extract_spatial","error, only mode 5 and can handle scalar input");
+ in = ggrd_open(argv[4],"r","hc_extract_spatial");
+ shps_read_d = hc_read_sh_solution(model,&dsol,in,binary_in,
+ verbose);
+ fclose(in);
+ ndsol = model->nradp2 * shps_read_d;
+ }
+ /*
+
+ room for spatial expansion
+
+ */
+ npoints = (vsol+i1*shps_read)->npoints;
+ if((vsol+i1*shps_read)->type != SH_RICK)
+ HC_ERROR("sh_extract_spatial","SH_RICK type required");
+ /* geographic set up */
+ nlat = (vsol+i1*shps_read)->rick.nlat;
+ nlon = (vsol+i1*shps_read)->rick.nlon;
+
+ ndata = npoints * shps ;
+ ndata_d = npoints * shps_read_d;
+ ndata_all = npoints * (shps + shps_read_d);
+
+ if((mode == 5)||(mode==6)){ /* save all layers */
+ hc_svecalloc(&data,model->nradp2 * ndata_all,"hc_extract_spatial");
+ }else
+ hc_svecalloc(&data, ndata_all,"hc_extract_spatial");
+ for(lc=0,ilayer=i1;ilayer <= i2;ilayer++,lc++){
/*
output
*/
@@ -111,30 +147,35 @@
if(verbose)
fprintf(stderr,"%s: printing v_r at layer %i (depth: %g)\n",argv[0],ilayer,zlabel);
- ivec=FALSE;sh_compute_spatial((sol+ilayer3),ivec,TRUE,&plm,data,verbose);
- sh_print_spatial_data_to_file((sol+ilayer3),shps,data,TRUE,zlabel,stdout);
+ ivec=FALSE;sh_compute_spatial((vsol+ilayer*shps_read),ivec,TRUE,&plm,data,verbose);
+ sh_print_spatial_data_to_file((vsol+ilayer*shps_read),shps,data,TRUE,zlabel,stdout);
break;
case 2:
/* */
if(verbose)
fprintf(stderr,"%s: printing v_theta v_phi SHE at layer %i (depth: %g)\n",argv[0],ilayer,zlabel);
- ivec=TRUE;sh_compute_spatial((sol+ilayer3+1),ivec,TRUE,&plm,data,verbose);
- sh_print_spatial_data_to_file((sol+ilayer3+1),shps,data,TRUE,zlabel,stdout);
+ ivec=TRUE;sh_compute_spatial((vsol+ilayer*shps_read+1),ivec,TRUE,&plm,data,verbose);
+ sh_print_spatial_data_to_file((vsol+ilayer*shps_read+1),shps,data,TRUE,zlabel,stdout);
break;
case 3:
if(verbose)
fprintf(stderr,"%s: printing v_r v_theta v_phi SHE at layer %i (depth: %g)\n",argv[0],ilayer,zlabel);
- ivec=FALSE;sh_compute_spatial((sol+ilayer3), ivec,TRUE,&plm,data,verbose); /* radial */
- ivec=TRUE; sh_compute_spatial((sol+ilayer3+1),ivec,TRUE,&plm,(data+npoints),verbose); /* theta,phi */
- sh_print_spatial_data_to_file((sol+ilayer3),shps,data,TRUE,zlabel,stdout);
+ ivec=FALSE;sh_compute_spatial((vsol+ilayer*shps_read), ivec,TRUE,&plm,data,verbose); /* radial */
+ ivec=TRUE; sh_compute_spatial((vsol+ilayer*shps_read+1),ivec,TRUE,&plm,(data+npoints),verbose); /* theta,phi */
+ sh_print_spatial_data_to_file((vsol+ilayer*shps_read),shps,data,TRUE,zlabel,stdout);
break;
case 4:
fprintf(stdout,"%5i %11g\n",ilayer,HC_Z_DEPTH(model->r[ilayer]));
break;
case 5: /* compute all and store */
case 6:
- ivec=FALSE;sh_compute_spatial((sol+ilayer3), ivec,TRUE,&plm,(data+lc*ndata),verbose); /* radial */
- ivec=TRUE; sh_compute_spatial((sol+ilayer3+1),ivec,TRUE,&plm,(data+lc*ndata+npoints),verbose); /* theta,phi */
+ ivec=FALSE;sh_compute_spatial((vsol+ilayer*shps_read), ivec,TRUE,&plm,(data+lc*ndata_all),verbose); /* radial */
+ ivec=TRUE; sh_compute_spatial((vsol+ilayer*shps_read+1),ivec,TRUE,&plm,(data+lc*ndata_all+npoints),verbose); /* theta,phi */
+ if(read_dsol){
+ if(!shps_read_d)
+ HC_ERROR("sh_extract_spatial","logic error");
+ ivec=FALSE;sh_compute_spatial((dsol+ilayer*shps_read_d),ivec,TRUE,&plm,(data+lc*ndata_all+npoints*shps),verbose); /* radial */
+ }
break;
default:
fprintf(stderr,"%s: error, mode %i undefined\n",argv[0],mode);
@@ -143,17 +184,24 @@
}
}
/* clear and exit */
- sh_free_expansion(sol,nsol);
+ sh_free_expansion(vsol,nvsol);
+ if(read_dsol)
+ sh_free_expansion(dsol,ndsol);
free(plm);
/* */
if((mode == 5)||(mode==6)){
+ /*
+ print the already stored properties
+ */
if(shps != 3)HC_ERROR("hc_extract_spatial","shps has to be 3 for mode 5 and 6");
/* convert */
- hc_svecalloc(&xpos,model->nradp2 * ndata,"hc_extract_spatial"); /* ndata = npoints * 3 */
+ hc_svecalloc(&xpos,model->nradp2 * ndata,"hc_extract_spatial");
hc_svecalloc(&xvec,model->nradp2 * ndata,"hc_extract_spatial");
+ if(read_dsol)
+ hc_svecalloc(&xscalar,model->nradp2 * ndata_d,"hc_extract_spatial");
for(i=0;i < npoints;i++){ /* loop through all points */
/* lon lat coordinates */
- sh_get_coordinates((sol+i1*3),i,&lon,&lat);
+ sh_get_coordinates((vsol+i1*3),i,&lon,&lat);
theta = LAT2THETA(lat);phi = LON2PHI(lon);
xtmp[0] = xtmp[1] = sin(theta);
xtmp[0] *= cos(phi); /* x */
@@ -165,21 +213,29 @@
/* this is the slow data storage loop but it avoids
recomputing the polar basis vector */
poff = ilayer * ndata + i*shps; /* point offset */
- for(j=0;j<3;j++){
+ for(j=0;j < 3;j++){
xpos[poff+j] = xtmp[j] * model->r[ilayer]; /* cartesian coordinates */
}
/* data are stored a bit weirdly, this makes for lots of
jumping around in memory ... */
- pvec[0] = data[ilayer*ndata+i];
- pvec[1] = data[ilayer*ndata+npoints+i];
- pvec[2] = data[ilayer*ndata+npoints*2+i];
+ pvec[0] = data[ilayer*ndata_all + i];
+ pvec[1] = data[ilayer*ndata_all + npoints +i];
+ pvec[2] = data[ilayer*ndata_all + npoints*2+i];
lonlatpv2cv_with_base(pvec,polar_base,(xvec+poff));
+ /* assign scalar fata if any */
+ for(j=0;j < shps_read_d;j++)
+ xscalar[j * model->nradp2 * ndata_d + ilayer * npoints + i] =
+ data[ilayer * ndata_all + npoints*(shps+j) + i];
}
}
free(data);
/* print in VTK format */
- hc_print_vtk(stdout,xpos,xvec,npoints,model->nradp2,(mode==6));
+ hc_print_vtk(stdout,xpos,xvec,npoints,model->nradp2,(mode==6),
+ shps_read_d,xscalar,nlon,nlat);
free(xvec);free(xpos);
+ if(shps_read_d)
+ free(xscalar);
+
}else{
free(data);
}
Modified: mc/1D/hc/trunk/hc_filenames.h
===================================================================
--- mc/1D/hc/trunk/hc_filenames.h 2011-07-17 13:59:04 UTC (rev 18780)
+++ mc/1D/hc/trunk/hc_filenames.h 2011-07-18 06:12:07 UTC (rev 18781)
@@ -23,5 +23,7 @@
#define HC_POLSOL_FILE "psol.dat"
+
+
#define HC_LAYER_OUT_FILE "vdepth.dat" /* depth [km] output */
#define HC_GEOID_FILE "geoid.ab" /* geoid output file */
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2011-07-17 13:59:04 UTC (rev 18780)
+++ mc/1D/hc/trunk/hc_init.c 2011-07-18 06:12:07 UTC (rev 18781)
@@ -853,8 +853,10 @@
break;
}
if(save_orig_danom && (!hc->dens_init)){
- /* make a copy of the original density anomaly before applying
- depth dependent scaling, only done once per run */
+ /*
+ make a copy of the original density anomaly before applying
+ depth dependent scaling, only done once per run
+ */
hc_get_blank_expansions(&hc->dens_anom_orig,hc->inho+1,hc->inho,
"hc_assign_density");
for(i=0;i<hc->inho;i++){
@@ -870,7 +872,9 @@
*/
for(i=0;i < hc->inho;i++){
- /* depth dependent factor? */
+ /*
+ depth dependent factor?
+ */
local_scale = hc_find_dens_scale(hc->rden[i],hc->dens_scale,dd_dens_scale,rdf,sdf,ndf);
sh_scale_expansion((hc->dens_anom+i),local_scale);
if(verbose >= 2){
@@ -951,8 +955,6 @@
free(dbot);free(dtop);
} /* end layer structure part */
-
-
hc->dens_init = TRUE;
}
/*
@@ -1335,13 +1337,17 @@
}
-void hc_select_pvel(HC_PREC time, struct pvels *pvel,struct sh_lms *p, hc_boolean verbose)
+void hc_select_pvel(HC_PREC time, struct pvels *pvel,
+ struct sh_lms *p, hc_boolean verbose)
{
if(pvel->n == 1){
if(verbose){
fprintf(stderr,"hc_select_pvel: only one plate velocity loadeed, disregarding time argument\n");
}
- /* only one plate velocity */
+ /*
+ only one plate velocity
+
+ */
sh_copy_lms(pvel->p,p);
sh_copy_lms((pvel->p+1),(p+1));
Modified: mc/1D/hc/trunk/hc_input.c
===================================================================
--- mc/1D/hc/trunk/hc_input.c 2011-07-17 13:59:04 UTC (rev 18780)
+++ mc/1D/hc/trunk/hc_input.c 2011-07-18 06:12:07 UTC (rev 18781)
@@ -4,11 +4,14 @@
routines dealing with input, pass sol initialized as NULL
+returns shps, the number of expansions. should be 3 for velocities and
+tractions, and 1 for density anomalies
$Id: hc_input.c,v 1.8 2004/12/20 05:18:12 becker Exp $
*/
-void hc_read_sh_solution(struct hcs *hc, struct sh_lms **sol, FILE *in,
+int hc_read_sh_solution(struct hcs *hc,
+ struct sh_lms **sol, FILE *in,
hc_boolean binary, hc_boolean verbose)
{
int nset,ilayer,shps,lmax,type,ivec,nsol,i,os,n;
@@ -20,13 +23,13 @@
*/
n = os = 0;
- while(sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer,&nset,
- &zlabel,&ivec,in,FALSE,binary,
- verbose)){
+ while(sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer,
+ &nset,&zlabel,&ivec,in,
+ FALSE,binary,verbose)){
hc->sh_type = type;
- if((shps != 3)||(ilayer != n)){
- fprintf(stderr,"hc_read_sh_solution: error: shps %i ilayer %i n %i\n",
- shps,ilayer,n);
+ if(ilayer != n){
+ fprintf(stderr,"hc_read_sh_solution: error: ilayer %i n %i\n",
+ ilayer,n);
exit(-1);
}
if(n == 0){ /* initialize */
@@ -49,17 +52,25 @@
*/
sh_read_coefficients_from_file((*sol+os),shps,-1,in,
binary,unity,verbose);
- if(verbose)
- fprintf(stderr,"hc_read_sh_solution: z: %8.3f |r|: %12.5e |pol|: %12.5e |tor|: %12.5e\n",
- HC_Z_DEPTH(hc->r[ilayer]),
- sqrt(sh_total_power((*sol+os))),sqrt(sh_total_power((*sol+os+1))),
- sqrt(sh_total_power((*sol+os+2))));
-
+ if(verbose){
+ if(shps == 3)
+ fprintf(stderr,"hc_read_sh_solution: z: %8.3f |r|: %12.5e |pol|: %12.5e |tor|: %12.5e\n",
+ HC_Z_DEPTH(hc->r[ilayer]),
+ sqrt(sh_total_power((*sol+os))),sqrt(sh_total_power((*sol+os+1))),
+ sqrt(sh_total_power((*sol+os+2))));
+ else
+ fprintf(stderr,"hc_read_sh_solution: z: %8.3f |scalar|: %12.5e\n",
+ HC_Z_DEPTH(hc->r[ilayer]),sqrt(sh_total_power((*sol+os))));
+ }
n++;
os += shps;
}
if(n != nset)
HC_ERROR("hc_read_sh_solution","read error");
- if(verbose)
+ if(verbose){
fprintf(stderr,"hc_read_sh_solution: read %i solution layer\n",nset);
+ if(shps != 3)
+ fprintf(stderr,"hc_read_sh_solution: WARNING: only scalar field read\n");
+ }
+ return shps;
}
Modified: mc/1D/hc/trunk/hc_misc.c
===================================================================
--- mc/1D/hc/trunk/hc_misc.c 2011-07-17 13:59:04 UTC (rev 18780)
+++ mc/1D/hc/trunk/hc_misc.c 2011-07-18 06:12:07 UTC (rev 18781)
@@ -26,6 +26,15 @@
if(! (*x))
HC_MEMERROR(message);
}
+/* integer vector allocation */
+void hc_ivecalloc(int **x,int n,char *message)
+{
+ *x = (int *)malloc(sizeof(int)*(size_t)n);
+ if(! (*x))
+ HC_MEMERROR(message);
+}
+
+
/* general floating point vector allocation */
void hc_vecalloc(HC_PREC **x,int n,char *message)
{
@@ -266,3 +275,25 @@
polar_base[7]= cp;
polar_base[8]= 0.0;
}
+/*
+
+ given a sorted vector y (y0<y1<...<yn-1) with n elements, return
+ the weights f1 for element i1 and weight f2 for element i2 such
+ that the interpolation is at y1
+
+ */
+void hc_linear_interpolate(HC_PREC *y, int n, HC_PREC y1,
+ int *i1, int *i2,
+ HC_PREC *f1, HC_PREC *f2)
+{
+ int n1;
+ n1 = n-1;
+ *i2 = 0;
+ while((*i2 < n1) && (y[*i2] < y1))
+ *i2 += 1;
+ if(*i2 == 0)
+ *i2 = 1;
+ *i1 = *i2 - 1;
+ *f2 = (y1 - y[*i1])/(y[*i2]-y[*i1]);
+ *f1 = 1.0 - *f2;
+}
Modified: mc/1D/hc/trunk/hc_output.c
===================================================================
--- mc/1D/hc/trunk/hc_output.c 2011-07-17 13:59:04 UTC (rev 18780)
+++ mc/1D/hc/trunk/hc_output.c 2011-07-18 06:12:07 UTC (rev 18781)
@@ -396,15 +396,25 @@
*/
void hc_print_vtk(FILE *out,float *xloc,float *xvec,
- int npoints,int nlay,
- hc_boolean binary)
+ int npoints_orig,int nlay,
+ hc_boolean binary,int shps_d,
+ float *xscalar,int nlon, int nlat)
{
- int i,ilay,ndata,poff;
+ int i,ilay,ndata,poff,j,ndata_d,nele_lay,nele_x,nele_y,
+ npe,npe1,ncon[12],k,nleft,nlon_m1,npoints,
+ nele_lay_reg,nele_lay_pole,tl,tr,nlay_m1,
+ nele_brick = 8,nele_tri = 6;
+
hc_boolean little_endian;
+ float xtmp[3],r,spole[3],npole[3];
/* determine machine type */
little_endian = hc_is_little_endian();
/* */
- ndata = npoints*3;
+ npoints = npoints_orig + 2;
+ ndata = npoints_orig*3;
+ ndata_d = npoints_orig * shps_d;
+
+ nlay_m1 = nlay - 1;
/* print VTK */
fprintf(out,"# vtk DataFile Version 4.0\n");
fprintf(out,"generated by hc_print_vtk\n");
@@ -412,11 +422,13 @@
fprintf(out,"BINARY\n");
else
fprintf(out,"ASCII\n");
- fprintf(out,"DATASET POLYDATA\n");
-
- fprintf(out,"POINTS %i float\n",npoints*nlay);
- for(ilay=0;ilay < nlay;ilay++){
- for(i=0;i < npoints;i++){
+ fprintf(out,"DATASET UNSTRUCTURED_GRID\n");
+ /*
+ nodes locations
+ */
+ fprintf(out,"POINTS %i float\n",npoints * nlay);
+ for(ilay=0;ilay < nlay;ilay++){ /* bottom up */
+ for(i=0;i < npoints_orig;i++){ /* S to N, W to E */
poff = ilay * ndata + i*3;
if(binary)
hc_print_be_float((xloc+poff),3,out,little_endian);
@@ -424,25 +436,209 @@
fprintf(out,"%.6e %.6e %.6e\n",
xloc[poff],xloc[poff+1],xloc[poff+2]);
}
+ /*
+ south and north poles, add two, to go to npoints per layer
+ */
+ poff = ilay * ndata;
+ r = sqrt(xloc[poff]*xloc[poff] +
+ xloc[poff+1]*xloc[poff+1] +
+ xloc[poff+2]*xloc[poff+2]);
+ /* south pole */
+ xtmp[0] = xtmp[1] = 0.0;xtmp[2]=-r;
+ if(binary)
+ hc_print_be_float((xtmp),3,out,little_endian);
+ else
+ fprintf(out,"%.6e %.6e %.6e\n",xtmp[0],xtmp[1],xtmp[2]);
+ /* north pole */
+ xtmp[2] = r;
+ if(binary)
+ hc_print_be_float((xtmp),3,out,little_endian);
+ else
+ fprintf(out,"%.6e %.6e %.6e\n",xtmp[0],xtmp[1],xtmp[2]);
}
+ /* */
+ nele_x = nlon;nlon_m1=nlon-1; /* elements in longitude */
+ nele_y = (nlat - 1); /* proper elements and pole connections */
+ /* top row node names */
+ tl = (nlat-1)*nlon;tr = tl + nlon;
+
+ nele_lay_reg = nele_x * nele_y; /* regular elements per layer */
+ nele_lay_pole = 2 * nele_x;
+ nele_lay = nele_lay_reg + nele_lay_pole; /* total per layer */
+ /*
+ element connectivity
+ */
+ fprintf(out,"CELLS %i %i\n",nlay_m1 * nele_lay,
+ nlay_m1 * (nele_lay_reg * (1+nele_brick)
+ + nele_x * 2 * (1+nele_tri)));
+ for(ilay = 0; ilay < nlay_m1; ilay++){
+ /*
+ loop bottom up, until one layer below top
+ */
+ npe = nele_brick; /* real nodes per element */
+ npe1 = npe+1;
+ ncon[0] = npe; /* counter */
+ for(i=0;i < nele_y;i++){
+ for(j=0;j < nele_x;j++){
+ nleft = ilay * npoints + i * nlon + j;
+ if(j == nlon_m1){ /* at edge, wrap around */
+ ncon[4] = nleft+nlon;ncon[3] = nleft-nlon_m1+nlon;
+ ncon[1] = nleft; ncon[2] = nleft-nlon_m1;
+ }else{ /* regular */
+ ncon[4] = nleft+nlon;ncon[3] = nleft+nlon+1;
+ ncon[1] = nleft; ncon[2] = nleft+1;
+ }
+ /* top level */
+ for(k=0;k < 4;k++)
+ ncon[5+k] = ncon[k+1] + npoints;
+ if(binary){
+ hc_print_be_int(ncon,npe1,out,little_endian);
+ }else{
+ for(k=0;k < npe1;k++)
+ fprintf(out,"%i ",ncon[k]);
+ fprintf(out,"\n");
+ }
+ }
+ }
+ /*
+ south and north polar rows
+ */
+ npe = nele_tri;npe1 = npe+1;
+ ncon[0] = npe;
+ for(j=0;j < 2;j++){
+ for(k=0;k < nele_x;k++){
+ if(j == 0){
+ /* south pole */
+ nleft = ilay * npoints + k;
+ ncon[1] = ilay * npoints + npoints_orig;
+ if(k == nlon_m1){ /* at edge, wrap around */
+ ncon[2] = nleft-nlon_m1;ncon[3]=nleft;
+ }else{
+ ncon[2] = nleft+1;ncon[3]=nleft;
+ }
+ }else{ /* north pole */
+ nleft = ilay * npoints + (nele_y-1) * nlon + k;
+ if(k == nlon_m1){ /* at edge, wrap around */
+ ncon[1] = nleft;ncon[2] = nleft-nlon_m1;
+ }else{
+ ncon[1] = nleft;ncon[2] = nleft+1;
+ }
+ ncon[3] = ilay * npoints + npoints_orig+1;
+ }
+ for(i=0;i < 3;i++)
+ ncon[4+i] = ncon[i+1] + npoints;
+ if(binary){
+ hc_print_be_int(ncon,npe1,out,little_endian);
+ }else{
+ for(i=0;i < npe1;i++)
+ fprintf(out,"%i ",ncon[i]);
+ fprintf(out,"\n");
+ }
+ }
+ }
+ }
+ /*
+ print cell types
+ */
+ fprintf(out,"CELL_TYPES %i\n",nlay_m1 * nele_lay);
+ for(ilay = 0; ilay < nlay_m1; ilay++){
+ if(binary){
+ /* binary */
+ ncon[0] = 12; /* VTK quad */
+ for(i=0;i < nele_lay_reg;i++)
+ hc_print_be_int(ncon,1,out,little_endian);
+ ncon[0] = 13; /* VTK triangle */
+ for(i=0;i < nele_lay_pole;i++)
+ hc_print_be_int(ncon,1,out,little_endian);
+ }else{
+ /* ascicc */
+ for(i=0;i < nele_lay_reg;i++){
+ fprintf(out,"%i ",12); /* vtk quad */
+ if(i % 80 == 0)
+ fprintf(out,"\n");
+ }
+ for(i=0;i < nele_lay_pole;i++){
+ fprintf(out,"%i ",13); /* vtk triagnle */
+ if(i % 80 == 0)
+ fprintf(out,"\n");
+ }
+ fprintf(out,"\n");
+ }
+ }
fprintf(out,"POINT_DATA %i\n",npoints*nlay);
+ if(shps_d){
+ for(j=0;j < shps_d;j++){
+ fprintf(out,"SCALARS scalar%i float 1\n",j+1);
+ fprintf(out,"LOOKUP_TABLE default\n");
+ for(ilay=0;ilay < nlay;ilay++){
+ spole[0] = npole[0] = 0.0; /* pole avg */
+ for(i=0,poff = j * nlay * ndata_d + ilay * ndata_d;
+ i < npoints_orig;i++,poff++){ /* all points */
+ if(binary)
+ hc_print_be_float((xscalar+poff),1,out,little_endian);
+ else{
+ fprintf(out,"%.6e ",xscalar[poff]);
+ if(i%20 == 0)fprintf(out,"\n");
+ }
+ if(i < nlon)
+ spole[0] += xscalar[poff];
+ if((i >= tl) && (i < tr))
+ npole[0] += xscalar[poff];
+ }
+ spole[0] /= (float)nlon;
+ npole[0] /= (float)nlon;
+ if(!binary){ /* ascii */
+ fprintf(out,"\n");
+ fprintf(out,"%.6e %.6e\n",spole[0],npole[0]);
+ }else{ /* binary */
+ hc_print_be_float(spole,1,out,little_endian);
+ hc_print_be_float(npole,1,out,little_endian);
+ }
+ }
+ }
+ }
fprintf(out,"VECTORS velocity float\n");
for(ilay=0;ilay < nlay;ilay++){
- for(i=0;i < npoints;i++){
+ spole[0] = spole[1] = spole[2] =
+ npole[0] = npole[1] = npole[2] = 0.0; /* pole avg */
+ for(i=0;i < npoints_orig;i++){
poff = ilay * ndata + i*3;
- if(binary)
+ if(i < nlon){
+ for(k=0;k<3;k++)
+ spole[k] += xvec[poff+k];
+ }
+ if((i >= tl) && (i < tr)){
+ for(k=0;k<3;k++)
+ npole[k] += xvec[poff+k];
+ }
+ if(binary) /* binary */
hc_print_be_float((xvec+poff),3,out,little_endian);
- else
+ else /* ascii */
fprintf(out,"%.6e %.6e %.6e\n",xvec[poff],xvec[poff+1],
xvec[poff+2]);
}
+ for(k=0;k<3;k++){
+ spole[k] /= (float)nlon;
+ npole[k] /= (float)nlon;
+ }
+ if(binary){
+ /* binary */
+ hc_print_be_float(spole,3,out,little_endian);
+ hc_print_be_float(npole,3,out,little_endian);
+ }else{
+ /* ascii */
+ fprintf(out,"%.6e %.6e %.6e\n",spole[0],spole[1],spole[2]);
+ fprintf(out,"%.6e %.6e %.6e\n",npole[0],npole[1],npole[2]);
+ }
}
}
/*
- print big endian binary
+ print big endian binary to file, no matter what hardware
+
*/
-void hc_print_be_float(float *x, int n, FILE *out, hc_boolean little_endian)
+void hc_print_be_float(float *x, int n, FILE *out,
+ hc_boolean little_endian)
{
int i;
float *xcopy;
@@ -460,6 +656,27 @@
fwrite(x,len,n,out);
}
}
+void hc_print_be_int(int *x, int n, FILE *out,
+ hc_boolean little_endian)
+{
+ int i, *xcopy;
+ static size_t len = sizeof(int);
+ if(little_endian){
+ /* need to flip the byte order */
+ hc_ivecalloc(&xcopy,n,"hc_print_be_int");
+ memcpy(xcopy,x,len*n);
+ for(i=0;i < n;i++)
+ hc_flip_byte_order((void *)(xcopy+i),len);
+ fwrite(xcopy,len,n,out);
+ free(xcopy);
+ }else{
+ /* can write as is */
+ fwrite(x,len,n,out);
+ }
+}
+/*
+ check if we're on a little endian machine
+ */
hc_boolean hc_is_little_endian(void)
{
static const unsigned long a = 1;
@@ -468,7 +685,7 @@
/*
-flip endianness of x
+ flip endianness of x
*/
void hc_flip_byte_order(void *x, size_t len)
@@ -483,7 +700,12 @@
hc_flipit(x,copy,len);
free(copy);
}
-/* this should not be called with (i,i,size i) */
+/*
+
+ actually flip the big endianness this should not be called with
+ (i,i,size i)
+
+*/
void hc_flipit(void *d, void *s, size_t len)
{
unsigned char *dest = d;
@@ -492,4 +714,32 @@
for (; len; len--)
*dest++ = *src--;
}
+/* print the density anomaly field interpolated to the nodal radii */
+void hc_print_dens_anom(struct hcs *hc,
+ FILE *out, hc_boolean binary,
+ hc_boolean verbose)
+{
+ int i,i1,i2;
+ HC_PREC f1,f2;
+ HC_PREC fac[3] = {1.,1.,1.};
+ struct sh_lms *exp;
+ sh_allocate_and_init(&exp,3,hc->dens_anom[0].lmax,hc->sh_type,0,FALSE,FALSE);
+ for(i=0;i < hc->nradp2;i++){
+ /* interpolate density depth to velocity node layer depth */
+ hc_linear_interpolate(hc->rden,hc->inho,hc->r[i],&i1,&i2,&f1,&f2);
+
+ sh_copy_lms((hc->dens_anom+i1),(exp+0));sh_scale_expansion((exp+0),f1);
+ sh_copy_lms((hc->dens_anom+i2),(exp+1));sh_scale_expansion((exp+1),f2);
+ sh_c_is_a_plus_b_coeff((exp+2),(exp+0),(exp+1)); /* c = a+b */
+ /* print to file */
+ sh_print_parameters_to_file((exp+2),1,i,hc->nradp2,HC_Z_DEPTH(hc->r[i]),out,FALSE,binary,verbose);
+ sh_print_coefficients_to_file((exp+2),1,out,fac,binary,verbose);
+ if(verbose>2)fprintf(stderr,"hc_print_dens_anom: z: %8.3f (f1: %6.3f f2: %6.3f) %3i/%3i pow: %10.3e %10.3e %10.3e\n",
+ HC_Z_DEPTH(hc->r[i]),f1,f2,i+1,hc->nradp2,
+ sqrt(sh_total_power((hc->dens_anom+i1))),sqrt(sh_total_power((hc->dens_anom+i2))),
+ sqrt(sh_total_power((exp+2))));
+ }
+
+ sh_free_expansion(exp,3);
+}
Modified: mc/1D/hc/trunk/main.c
===================================================================
--- mc/1D/hc/trunk/main.c 2011-07-17 13:59:04 UTC (rev 18780)
+++ mc/1D/hc/trunk/main.c 2011-07-18 06:12:07 UTC (rev 18781)
@@ -194,6 +194,24 @@
p->solution_mode,
p->sol_binary_out,p->verbose);
fclose(out);
+ /*
+ print the density field
+ */
+ sprintf(file_prefix,"dscaled");
+ if(p->sol_binary_out)
+ sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_BINARY);
+ else
+ sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_ASCII);
+ if(p->verbose)
+ fprintf(stderr,"%s: writing scaled density anomaly field to %s\n",
+ argv[0],filename);
+
+ out = ggrd_open(filename,"w","main");
+ hc_print_dens_anom(model,out,p->sol_binary_out,p->verbose);
+ fclose(out);
+
+
+ /* */
if(p->compute_geoid){
if(p->compute_geoid_correlations){
if(p->compute_geoid == 2){ /* check if all geoids were computed */
Modified: mc/1D/hc/trunk/sh.h
===================================================================
--- mc/1D/hc/trunk/sh.h 2011-07-17 13:59:04 UTC (rev 18780)
+++ mc/1D/hc/trunk/sh.h 2011-07-18 06:12:07 UTC (rev 18781)
@@ -76,7 +76,8 @@
*/
hc_boolean plm_computed,plm_computed_irr;
- int old_lmax,old_ivec,old_tnplm,old_tnplm_irr,old_lmax_irr,old_ivec_irr;
+ int old_lmax,old_ivec,old_tnplm,old_tnplm_irr,
+ old_lmax_irr,old_ivec_irr;
/*
Modified: mc/1D/hc/trunk/sh_exp.c
===================================================================
--- mc/1D/hc/trunk/sh_exp.c 2011-07-17 13:59:04 UTC (rev 18780)
+++ mc/1D/hc/trunk/sh_exp.c 2011-07-18 06:12:07 UTC (rev 18781)
@@ -1955,13 +1955,16 @@
}
}
-/* copy a whole expansion structure */
+/*
+ copy a whole expansion structure
+ a --> b, i.e.
+ b = a
+
+
+*/
void sh_copy_lms(struct sh_lms *a, struct sh_lms *b)
{
- memcpy(b,a,sizeof(struct sh_lms));
-
- /* make sure this really gets done, if memory is not contiguous */
b->type = a->type;
b->lmax = a->lmax;
b->spectral_init = a->spectral_init;
@@ -1981,8 +1984,12 @@
b->old_tnplm_irr = a->old_tnplm_irr;
b->old_lmax_irr = a->old_lmax_irr;
b->old_ivec_irr = a->old_ivec_irr;
-
- sh_aexp_equals_bexp_coeff(a,b);
+ /* this still needs to be fixed */
+ b->rick.nlat = a->rick.nlat;
+ b->rick.nlon = a->rick.nlon;
+
+ sh_aexp_equals_bexp_coeff(b,a);
+
}
/*
copy the coefficients of one expansion to another
@@ -2036,6 +2043,49 @@
break;
}
}
+/* sum two expansions, c = a + b */
+
+void sh_c_is_a_plus_b_coeff(struct sh_lms *c, struct sh_lms *a, struct sh_lms *b)
+{
+ int i;
+ if((a->type != b->type)||(b->type != c->type)){
+ fprintf(stderr,"sh_c_is_a_plus_b_coeff: error: type mix (%i vs. %i vs. %i) not implemented yet\n",
+ a->type,b->type,c->type);
+ exit(-1);
+ }
+ /*
+ make sure lmax(b) = lmax(a) = lmax(c)
+
+ */
+ if((a->lmax != b->lmax)||(b->lmax != c->lmax)){
+ fprintf(stderr,"sh_c_is_a_plus_b_coeff: error: a!=b || b != c lmax %i %i %i\n",
+ a->lmax,b->lmax,c->lmax);
+ exit(-1);
+ }
+ switch(a->type){
+#ifdef HC_USE_HEALPIX
+
+ case SH_HEALPIX: /* nplm should reflect the lmax */
+ for(i=0;i < b->n_lm;i++){
+ c->alm_c[i].dr = a->alm_c[i].dr + b->alm_c[i].dr;
+ c->alm_c[i].dr = a->alm_c[i].di + b->alm_c[i].di;
+ }
+ break;
+#endif
+ case SH_RICK:
+ for(i=0;i < b->n_lm;i++)
+ c->alm[i] = a->alm[i] + b->alm[i];
+ break;
+#ifdef HC_USE_SPHEREPACK
+ case SH_SPHEREPACK_GAUSS:
+ case SH_SPHEREPACK_EVEN:
+ break;
+#endif
+ default:
+ sh_exp_type_error("sh_c_is_a_plus_b_coeff",a);
+ break;
+ }
+}
/*
given a spherical harmonic expansion, scale it with factor fac[0...lmax]
More information about the CIG-COMMITS
mailing list