[cig-commits] r20528 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Wed Jul 18 15:19:21 PDT 2012
Author: becker
Date: 2012-07-18 15:19:20 -0700 (Wed, 18 Jul 2012)
New Revision: 20528
Modified:
mc/1D/hc/trunk/Makefile
mc/1D/hc/trunk/Makefile.include
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/hc_polsol.c
mc/1D/hc/trunk/make_tar
Log:
Implemented approximation for CMB mechanical boundary condition
and location for high degrees to improve stability following
Steinberger & Torsvik
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2012-07-16 12:55:19 UTC (rev 20527)
+++ mc/1D/hc/trunk/Makefile 2012-07-18 22:19:20 UTC (rev 20528)
@@ -258,7 +258,7 @@
mkdir -p $(BDIR);
clean:
- rm -f $(ODIR)/*.o $(ODIR)/*.a $(BDIR)/*
+ rm -f hc_auto_proto.h $(ODIR)/*.o $(ODIR)/*.a $(BDIR)/*
#
# library
Modified: mc/1D/hc/trunk/Makefile.include
===================================================================
--- mc/1D/hc/trunk/Makefile.include 2012-07-16 12:55:19 UTC (rev 20527)
+++ mc/1D/hc/trunk/Makefile.include 2012-07-18 22:19:20 UTC (rev 20528)
@@ -6,6 +6,9 @@
#
# for GMT version >= 4.1.2, uncomment the next two lines
#
-#ADD_FLAGS = -Werror
+# quad precision
+ADD_FLAGS = -DHC_PRECISION=32 -O2
+#ADD_FLAGS = -O2
+
GGRD_INC_FLAGS = -I$(GMTHOME)/include -I$(NETCDFHOME)/include
GGRD_LIBS_LINKLINE = -lggrd -lgmt -lpsl -lnetcdf
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2012-07-16 12:55:19 UTC (rev 20527)
+++ mc/1D/hc/trunk/hc.h 2012-07-18 22:19:20 UTC (rev 20528)
@@ -175,6 +175,8 @@
abg_init , /* alpha, beta factors */
prop_mats_init; /* will be true only if save_prop_mats is
requested */
+ int solver_kludge_l; /* for treating higher orders
+ differently */
/* for solve */
hc_boolean tor_init, pol_init;
};
@@ -210,6 +212,8 @@
int pvel_mode; /* plate velocity mode */
HC_PREC pvel_time; /* time to use */
+ int solver_kludge_l; /* for CMB BC tricks */
+
hc_boolean solver_mode;
hc_boolean visc_init_mode;
HC_PREC elayer[4];
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2012-07-16 12:55:19 UTC (rev 20527)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2012-07-18 22:19:20 UTC (rev 20528)
@@ -40,86 +40,86 @@
void hc_struc_init(struct hcs **);
void hc_init_polsol_struct(struct hc_ps *);
void hc_init_main(struct hcs *, int, struct hc_parameters *);
-void hc_init_constants(struct hcs *, double, char *, unsigned short);
+void hc_init_constants(struct hcs *, long double, char *, unsigned short);
void hc_handle_command_line(int, char **, struct hc_parameters *);
-void hc_assign_viscosity(struct hcs *, int, double [4], struct hc_parameters *);
-void hc_assign_density(struct hcs *, unsigned short, int, char *, int, unsigned short, unsigned short, unsigned short, unsigned short, unsigned short, unsigned short, int, double *, double *, unsigned short);
-double hc_find_dens_scale(double, double, unsigned short, double *, double *, int);
+void hc_assign_viscosity(struct hcs *, int, long double [4], struct hc_parameters *);
+void hc_assign_density(struct hcs *, unsigned short, int, char *, int, unsigned short, unsigned short, unsigned short, unsigned short, unsigned short, unsigned short, int, long double *, long double *, unsigned short);
+long double hc_find_dens_scale(long double, long double, unsigned short, long double *, long double *, int);
void hc_init_phase_boundaries(struct hcs *, int, unsigned short);
void hc_assign_plate_velocities(struct hcs *, int, char *, unsigned short, int, unsigned short, unsigned short, unsigned short);
void hc_init_single_plate_exp(char *, struct hcs *, unsigned short, struct sh_lms *, unsigned short, unsigned short, unsigned short);
void hc_init_l_factors(struct hcs *, int);
void hc_get_blank_expansions(struct sh_lms **, int, int, char *);
void hc_struc_free(struct hcs **);
-void hc_assign_dd_scaling(int, double [4], struct hc_parameters *, double);
+void hc_assign_dd_scaling(int, long double [4], struct hc_parameters *, long double);
void hc_read_geoid(struct hc_parameters *);
-void hc_select_pvel(double, struct pvels *, struct sh_lms *, unsigned short);
+void hc_select_pvel(long double, struct pvels *, struct sh_lms *, unsigned short);
/* hc_input.c */
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 *);
+void hc_ludcmp_3x3(long double [3][3], int, int *);
+void hc_lubksb_3x3(long double [3][3], int, int *, long double *);
/* hc_misc.c */
-void hc_hvecalloc(double **, int, char *);
+void hc_hvecalloc(long double **, int, char *);
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_vecalloc(long double **, int, char *);
void hc_scmplx_vecalloc(struct hc_scmplx **, int, char *);
void hc_svecrealloc(float **, int, char *);
-void hc_dvecrealloc(double **, int, char *);
-void hc_vecrealloc(double **, int, char *);
-float hc_vec_rms_diff(double *, double *, int);
-float hc_vec_rms(double *, int);
+void hc_dvecrealloc(long double **, int, char *);
+void hc_vecrealloc(long double **, int, char *);
+float hc_vec_rms_diff(long double *, long double *, int);
+float hc_vec_rms(long double *, int);
void hc_a_equals_b_svector(float *, float *, int);
-void hc_a_equals_b_vector(double *, double *, int);
+void hc_a_equals_b_vector(long double *, long double *, int);
float hc_mean_svec(float *, int);
-double hc_mean_vec(double *, int);
-void hc_zero_dvector(double *, int);
+long double hc_mean_vec(long double *, int);
+void hc_zero_dvector(long double *, int);
void hc_zero_lvector(unsigned short *, int);
void hc_get_flt_frmt_string(char *, int, unsigned short);
char *hc_name_boolean(unsigned short);
unsigned short hc_toggle_boolean(unsigned short *);
void hc_advance_argument(int *, int, char **);
-void hc_compute_correlation(struct sh_lms *, struct sh_lms *, double *, int, unsigned short);
-void lonlatpv2cv(double, float, double *, double *);
-void lonlatpv2cv_with_base(double *, double *, double *);
-void calc_polar_base_at_theta_phi(double, double, double *);
-void hc_linear_interpolate(double *, int, double, int *, int *, double *, double *);
+void hc_compute_correlation(struct sh_lms *, struct sh_lms *, long double *, int, unsigned short);
+void lonlatpv2cv(long double, float, long double *, long double *);
+void lonlatpv2cv_with_base(long double *, long double *, long double *);
+void calc_polar_base_at_theta_phi(long double, long double, long double *);
+void hc_linear_interpolate(long double *, int, long double, int *, int *, long double *, long 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);
-void hc_print_spatial_solution(struct hcs *, struct sh_lms *, double *, char *, char *, int, unsigned short, unsigned short);
+void hc_print_spatial_solution(struct hcs *, struct sh_lms *, long double *, char *, char *, int, unsigned short, unsigned short);
void hc_print_depth_layers(struct hcs *, FILE *, unsigned short);
-void hc_print_3x3(double [3][3], FILE *);
-void hc_print_sm(double [6][4], FILE *);
-void hc_print_vector(double *, int, FILE *);
-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, double, double *);
+void hc_print_3x3(long double [3][3], FILE *);
+void hc_print_sm(long double [6][4], FILE *);
+void hc_print_vector(long double *, int, FILE *);
+void hc_print_vector_label(long double *, int, FILE *, char *);
+void hc_print_matrix_label(long double *, int, int, FILE *, char *);
+void hc_print_vector_row(long double *, int, FILE *);
+void hc_compute_solution_scaling_factors(struct hcs *, int, long double, long double, long 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 *, double *, double *, int, int, unsigned short, int, double *, int, int);
-int hc_print_be_float(double *, int, FILE *, unsigned short);
-int hc_print_float(double *, int, FILE *);
-int hc_read_float(double *, int, FILE *);
+void hc_print_toroidal_solution(long double *, int, struct hcs *, int, char *, unsigned short);
+void hc_print_vtk(FILE *, long double *, long double *, int, int, unsigned short, int, long double *, int, int);
+int hc_print_be_float(long double *, int, FILE *, unsigned short);
+int hc_print_float(long double *, int, FILE *);
+int hc_read_float(long double *, int, FILE *);
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);
+void hc_polsol(struct hcs *, int, long double *, int, long double *, unsigned short, struct sh_lms *, unsigned short, int, long double *, long double *, unsigned short, struct sh_lms *, struct sh_lms *, unsigned short, struct sh_lms *, unsigned short, unsigned short);
/* hc_propagator.c */
-void hc_evalpa(int, double, double, double, double *);
-void hc_evppot(int, double, double *);
+void hc_evalpa(int, long double, long double, long double, long double *);
+void hc_evppot(int, long double, long double *);
/* hc_solve.c */
void hc_solve(struct hcs *, unsigned short, int, struct sh_lms *, unsigned short, unsigned short, unsigned short, unsigned short, unsigned short, struct sh_lms *, struct sh_lms *, struct sh_lms *, unsigned short);
void hc_sum(struct hcs *, int, struct sh_lms *, struct sh_lms *, int, unsigned short, struct sh_lms *, unsigned short);
-void hc_compute_sol_spatial(struct hcs *, struct sh_lms *, double **, unsigned short);
+void hc_compute_sol_spatial(struct hcs *, struct sh_lms *, long double **, unsigned short);
/* hc_torsol.c */
-void hc_torsol(struct hcs *, int, int, int, double *, double **, double **, struct sh_lms *, struct sh_lms *, double *, unsigned short);
+void hc_torsol(struct hcs *, int, int, int, long double *, long double **, long double **, struct sh_lms *, struct sh_lms *, long double *, unsigned short);
/* main.c */
/* prem2dsm.c */
/* prem_util.c */
@@ -135,87 +135,87 @@
int prem_read_para_set(double *, int, int, FILE *);
/* print_gauss_lat.c */
/* rick_fft_c.c */
-void rick_cs2ab(double *, int);
-void rick_ab2cs(double *, int);
-void rick_realft_nr(double *, int, int);
-void rick_four1_nr(double *, int, int);
+void rick_cs2ab(long double *, int);
+void rick_ab2cs(long double *, int);
+void rick_realft_nr(long double *, int, int);
+void rick_four1_nr(long double *, int, int);
/* rick_sh_c.c */
-void rick_compute_allplm(int, int, double *, double *, struct rick_module *);
-void rick_compute_allplm_reg(int, int, double *, double *, struct rick_module *, double *, int);
-void rick_pix2ang(int, int, double *, double *, struct rick_module *);
-void rick_shc2d(double *, double *, int, int, double *, double *, struct rick_module *);
-void rick_shc2d_reg(double *, double *, int, int, double *, double *, struct rick_module *, double *, int, double *, int, unsigned short);
-void rick_shc2d_pre(double *, double *, int, double *, double *, int, double *, double *, struct rick_module *);
-void rick_shc2d_pre_reg(double *, double *, int, double *, double *, int, double *, double *, struct rick_module *, double *, int, double *, int, unsigned short);
-void rick_shc2d_irreg(double *, double *, int, int, double *, double *, struct rick_module *, double *, double *, int);
-void rick_shd2c(double *, double *, int, int, double *, double *, struct rick_module *);
-void rick_shd2c_pre(double *, double *, int, double *, double *, int, double *, double *, struct rick_module *);
+void rick_compute_allplm(int, int, long double *, long double *, struct rick_module *);
+void rick_compute_allplm_reg(int, int, long double *, long double *, struct rick_module *, long double *, int);
+void rick_pix2ang(int, int, long double *, long double *, struct rick_module *);
+void rick_shc2d(long double *, long double *, int, int, long double *, long double *, struct rick_module *);
+void rick_shc2d_reg(long double *, long double *, int, int, long double *, long double *, struct rick_module *, long double *, int, long double *, int, unsigned short);
+void rick_shc2d_pre(long double *, long double *, int, long double *, long double *, int, long double *, long double *, struct rick_module *);
+void rick_shc2d_pre_reg(long double *, long double *, int, long double *, long double *, int, long double *, long double *, struct rick_module *, long double *, int, long double *, int, unsigned short);
+void rick_shc2d_irreg(long double *, long double *, int, int, long double *, long double *, struct rick_module *, long double *, long double *, int);
+void rick_shd2c(long double *, long double *, int, int, long double *, long double *, struct rick_module *);
+void rick_shd2c_pre(long double *, long double *, int, long double *, long double *, int, long double *, long double *, struct rick_module *);
void rick_init(int, int, int *, int *, int *, struct rick_module *, unsigned short);
void rick_free_module(struct rick_module *, int);
-void rick_plmbar1(double *, double *, int, int, double, struct rick_module *);
-void rick_gauleg(double, double, double *, double *, int);
+void rick_plmbar1(long double *, long double *, int, int, long double, struct rick_module *);
+void rick_gauleg(long double, long double, long double *, long double *, int);
/* rotvec2vel.c */
FILE *rv_myopen(const char *, const char *);
/* sh_ana.c */
/* shana_sh.c */
void shana_compute_allplm(int, int, double *, double *, struct shana_module *);
void shana_pix2ang(int, int, double *, double *, struct shana_module *);
-void shana_shc2d(double *, double *, int, int, double *, double *, struct shana_module *);
-void shana_shc2d_pre(double *, double *, int, double *, double *, int, float *, float *, struct shana_module *);
-void shana_shd2c(double *, double *, int, int, double *, double *, struct shana_module *);
-void shana_shd2c_pre(double *, double *, int, double *, double *, int, double *, double *, struct shana_module *);
+void shana_shc2d(long double *, long double *, int, int, long double *, long double *, struct shana_module *);
+void shana_shc2d_pre(long double *, long double *, int, double *, double *, int, float *, float *, struct shana_module *);
+void shana_shd2c(long double *, long double *, int, int, long double *, long double *, struct shana_module *);
+void shana_shd2c_pre(long double *, long double *, int, double *, double *, int, long double *, long double *, struct shana_module *);
void shana_init(int, int, int *, int *, int *, struct shana_module *);
void shana_free_module(struct shana_module *, int);
-void shana_plmbar1(double *, double *, int, int, double, struct shana_module *);
+void shana_plmbar1(double *, double *, int, int, long double, struct shana_module *);
/* sh_corr.c */
/* sh_exp.c */
void sh_allocate_and_init(struct sh_lms **, int, int, int, int, unsigned short, unsigned short);
void sh_init_expansion(struct sh_lms *, int, int, int, unsigned short, unsigned short);
void sh_free_expansion(struct sh_lms *, int);
void sh_clear_alm(struct sh_lms *);
-double sh_total_power(struct sh_lms *);
-void sh_compute_power_per_degree(struct sh_lms *, double *);
-double sh_correlation(struct sh_lms *, struct sh_lms *, int);
-double sh_correlation_per_degree(struct sh_lms *, struct sh_lms *, int, int);
-void sh_print_parameters_to_file(struct sh_lms *, int, int, int, double, FILE *, unsigned short, unsigned short, unsigned short);
-unsigned short sh_read_parameters_from_file(int *, int *, int *, int *, int *, double *, int *, FILE *, unsigned short, unsigned short, unsigned short);
-void sh_print_coefficients_to_file(struct sh_lms *, int, FILE *, double *, unsigned short, unsigned short);
-void sh_read_coefficients_from_file(struct sh_lms *, int, int, FILE *, unsigned short, double *, unsigned short);
+long double sh_total_power(struct sh_lms *);
+void sh_compute_power_per_degree(struct sh_lms *, long double *);
+long double sh_correlation(struct sh_lms *, struct sh_lms *, int);
+long double sh_correlation_per_degree(struct sh_lms *, struct sh_lms *, int, int);
+void sh_print_parameters_to_file(struct sh_lms *, int, int, int, long double, FILE *, unsigned short, unsigned short, unsigned short);
+unsigned short sh_read_parameters_from_file(int *, int *, int *, int *, int *, long double *, int *, FILE *, unsigned short, unsigned short, unsigned short);
+void sh_print_coefficients_to_file(struct sh_lms *, int, FILE *, long double *, unsigned short, unsigned short);
+void sh_read_coefficients_from_file(struct sh_lms *, int, int, FILE *, unsigned short, long double *, unsigned short);
void sh_print_nonzero_coeff(struct sh_lms *, FILE *);
-void sh_read_spatial_data_from_file(struct sh_lms *, FILE *, unsigned short, int, double *, double *);
-void sh_read_spatial_data_from_grd(struct sh_lms *, struct ggrd_gt *, unsigned short, int, double *, double *);
-void sh_read_spatial_data(struct sh_lms *, FILE *, struct ggrd_gt *, unsigned short, unsigned short, int, double *, double *);
-void sh_compute_spatial_basis(struct sh_lms *, FILE *, unsigned short, double, double **, int, unsigned short);
-void sh_compute_spectral(double *, int, unsigned short, double **, struct sh_lms *, unsigned short);
-void sh_compute_spatial(struct sh_lms *, int, unsigned short, double **, double *, unsigned short);
-void sh_compute_spatial_reg(struct sh_lms *, int, unsigned short, double **, double *, int, double *, int, double *, unsigned short, unsigned short);
-void sh_compute_spatial_irreg(struct sh_lms *, int, double *, double *, int, double *, unsigned short);
+void sh_read_spatial_data_from_file(struct sh_lms *, FILE *, unsigned short, int, long double *, long double *);
+void sh_read_spatial_data_from_grd(struct sh_lms *, struct ggrd_gt *, unsigned short, int, long double *, long double *);
+void sh_read_spatial_data(struct sh_lms *, FILE *, struct ggrd_gt *, unsigned short, unsigned short, int, long double *, long double *);
+void sh_compute_spatial_basis(struct sh_lms *, FILE *, unsigned short, long double, long double **, int, unsigned short);
+void sh_compute_spectral(long double *, int, unsigned short, long double **, struct sh_lms *, unsigned short);
+void sh_compute_spatial(struct sh_lms *, int, unsigned short, long double **, long double *, unsigned short);
+void sh_compute_spatial_reg(struct sh_lms *, int, unsigned short, long double **, long double *, int, long double *, int, long double *, unsigned short, unsigned short);
+void sh_compute_spatial_irreg(struct sh_lms *, int, long double *, long double *, int, long double *, unsigned short);
void sh_exp_type_error(char *, struct sh_lms *);
-void sh_print_plm(double *, int, int, int, FILE *);
-void sh_print_spatial_data_to_file(struct sh_lms *, int, double *, unsigned short, double, FILE *);
-void sh_get_coordinates(struct sh_lms *, int, double *, double *);
-void sh_print_reg_spatial_data_to_file(struct sh_lms *, int, double *, unsigned short, double, double *, int, double *, int, FILE *);
-void sh_print_irreg_spatial_data_to_file(struct sh_lms *, int, double *, unsigned short, double, double *, double *, int, FILE *);
-void sh_compute_plm(struct sh_lms *, int, double **, unsigned short);
-void sh_compute_plm_reg(struct sh_lms *, int, double **, unsigned short, double *, int);
-void sh_get_coeff(struct sh_lms *, int, int, int, unsigned short, double *);
-void sh_write_coeff(struct sh_lms *, int, int, int, unsigned short, double *);
-void sh_add_coeff(struct sh_lms *, int, int, int, unsigned short, double *);
+void sh_print_plm(long double *, int, int, int, FILE *);
+void sh_print_spatial_data_to_file(struct sh_lms *, int, long double *, unsigned short, long double, FILE *);
+void sh_get_coordinates(struct sh_lms *, int, long double *, long double *);
+void sh_print_reg_spatial_data_to_file(struct sh_lms *, int, long double *, unsigned short, long double, long double *, int, long double *, int, FILE *);
+void sh_print_irreg_spatial_data_to_file(struct sh_lms *, int, long double *, unsigned short, long double, long double *, long double *, int, FILE *);
+void sh_compute_plm(struct sh_lms *, int, long double **, unsigned short);
+void sh_compute_plm_reg(struct sh_lms *, int, long double **, unsigned short, long double *, int);
+void sh_get_coeff(struct sh_lms *, int, int, int, unsigned short, long double *);
+void sh_write_coeff(struct sh_lms *, int, int, int, unsigned short, long double *);
+void sh_add_coeff(struct sh_lms *, int, int, int, unsigned short, long 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);
+void sh_scale_expansion_l_factor(struct sh_lms *, long double *);
+void sh_scale_expansion(struct sh_lms *, long double);
/* sh_extract_layer.c */
/* sh_model.c */
void sh_init_model(struct sh_lms_model *, int, int, int, int, int, int, unsigned short);
void sh_free_model(struct sh_lms_model *);
void sh_print_model_coefficients(struct sh_lms_model *, FILE *, unsigned short, unsigned short);
void sh_print_model_spatial_basis(struct sh_lms_model *, FILE *, unsigned short);
-void sh_read_model_spatial_data(struct sh_lms_model *, double **, FILE *, unsigned short);
-void sh_compute_model_spectral(struct sh_lms_model *, double *, unsigned short);
-void sh_compute_model_spatial(struct sh_lms_model *, double **, unsigned short);
-void sh_print_model_spatial_data(struct sh_lms_model *, double *, FILE *, unsigned short);
+void sh_read_model_spatial_data(struct sh_lms_model *, long double **, FILE *, unsigned short);
+void sh_compute_model_spectral(struct sh_lms_model *, long double *, unsigned short);
+void sh_compute_model_spatial(struct sh_lms_model *, long double **, unsigned short);
+void sh_print_model_spatial_data(struct sh_lms_model *, long double *, FILE *, unsigned short);
/* sh_power.c */
/* sh_syn.c */
/* sh_test.c */
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2012-07-16 12:55:19 UTC (rev 20527)
+++ mc/1D/hc/trunk/hc_init.c 2012-07-18 22:19:20 UTC (rev 20528)
@@ -60,7 +60,7 @@
p->elayer[0] = 50.; p->elayer[1] = 1.; p->elayer[2] = 0.1; p->elayer[3] = 50.;
p->visc_init_mode = HC_INIT_E_FROM_FILE; /* by default, read viscosity from file */
-
+ p->solver_kludge_l = INT_MAX; /* default: no solver tricks */
/*
depth dependent scaling of density files?
*/
@@ -117,6 +117,9 @@
psp->abg_init = FALSE; /* alpha, beta factors */
psp->prop_mats_init = FALSE; /* will be true only if save_prop_mats is */
psp->tor_init = psp->pol_init = FALSE;
+ psp->solver_kludge_l = INT_MAX; /* every l > psp->solver_kludge_l will
+ have modified core boundary
+ conditions */
}
/*
@@ -243,7 +246,12 @@
}else{
HC_ERROR("hc_init","boundary condition logic error");
}
-
+ /* solver tricks */
+ hc->psp.solver_kludge_l = p->solver_kludge_l;
+ if(hc->psp.solver_kludge_l != INT_MAX)
+ if(p->verbose)
+ fprintf(stderr,"hc_init_main: WARNING: applying solver CMB kludge for l > %i\n",hc->psp.solver_kludge_l);
+
/*
phase boundaries, if any
*/
@@ -413,6 +421,8 @@
(double)p->pvel_time);
fprintf(stderr,"solution procedure and I/O options:\n");
+ fprintf(stderr,"-cbckl\tval\twill modify CMB boundary condition for all l > val with solver kludge (%i)\n",
+ p->solver_kludge_l);
fprintf(stderr,"-ng\t\tdo not compute and print the geoid (%i)\n",
p->compute_geoid);
fprintf(stderr,"-ag\t\tcompute geoid at all layer depths, as opposed to the surface only\n");
@@ -432,6 +442,9 @@
}else if(strcmp(argv[i],"-ds")==0){ /* density anomaly scaling factor */
hc_advance_argument(&i,argc,argv);
sscanf(argv[i],HC_FLT_FORMAT,&p->dens_anom_scale);
+ }else if(strcmp(argv[i],"-cbckl")==0){ /* solver kludge */
+ hc_advance_argument(&i,argc,argv);
+ sscanf(argv[i],"%i",&p->solver_kludge_l);
}else if(strcmp(argv[i],"-vtime")==0){ /* */
hc_advance_argument(&i,argc,argv);
sscanf(argv[i],HC_FLT_FORMAT,&p->pvel_time);
Modified: mc/1D/hc/trunk/hc_polsol.c
===================================================================
--- mc/1D/hc/trunk/hc_polsol.c 2012-07-16 12:55:19 UTC (rev 20527)
+++ mc/1D/hc/trunk/hc_polsol.c 2012-07-18 22:19:20 UTC (rev 20528)
@@ -196,13 +196,14 @@
nprops_max,jsol;
double *xprem;
HC_HIGH_PREC *b,du1,du2,el,rnext,drho,dadd;
+ HC_PREC rbound_kludge;
HC_HIGH_PREC amat[3][3],bvec[3],u[4],poten[2],
unew[4],potnew[2],clm[2];
/*
structures which hold u[6][4] type arrays
*/
struct hc_sm cmb, *u3;
- hc_boolean qvis,qinho,hit;
+ hc_boolean qvis,qinho,hit,kludge_warned;
/*
define a few offset and size pointers
*/
@@ -524,6 +525,19 @@
*/
el = (HC_PREC)l;
+
+ /*
+ this will normally be a very small number so that all
+ propagators will be computed above the regular CMB
+
+ if solver_kludge_l is set to within the [0;L] domain, the depth
+ of the bottom will depend on l
+
+ (rprops(i).ge.(1.-(1.-0.5448)*50./l))
+ */
+ rbound_kludge = (1. - (1.-hc->r_cmb)*(HC_PREC)hc->psp.solver_kludge_l/el);
+ kludge_warned = FALSE;
+
if((!save_prop_mats) || (!hc->psp.prop_mats_init)|| (viscosity_or_layer_changed)){
//
// get all propagators now, as they only depend on l
@@ -592,12 +606,25 @@
for(ibv=0;ibv < 4;ibv++)
cmb.u[i6][ibv] = 0.0;
- cmb.u[1][1] = 1.0; /* set this to zero for no-slip,
+ if(l > hc->psp.solver_kludge_l){
+ /*
+ solver trick to ensure stabilty following Steinberger &
+ Torsvik, doi:10.1029/2011GC003808
+
+ ucmb(4,1)=1.d0
+ */
+ cmb.u[3][1] = 1.0; /* make core fixed */
+ }else{
+ /* regular operation */
+ /* ucmb(2,1)=1.d0 */
+ cmb.u[1][1] = 1.0; /* set this to zero for no-slip,
in general CMB is free slip */
- cmb.u[2][2] = 1.0;
- cmb.u[4][3] = 1.0;
- cmb.u[5][3] = el;
+ }
+ cmb.u[2][2] = 1.0; /* ucmb(3,2)=1.d0 */
+ cmb.u[4][3] = 1.0; /* ucmb(5,3)=1.d0 */
+ cmb.u[5][3] = el; /* ucmb(6,3)=float(l) */
+
for(ibv=0;ibv < 4;ibv++){
/*
@@ -637,67 +664,79 @@
I NPROPS LOOP
*/
+
- //
- //
- // PROPAGATE U TO NEXT RADIUS IN RPROPS
- //
- for(os=pos1 + i*16,i2=0;i2 < 4;i2++,os += 4){
- unew[i2] = 0.0;
- for(i3=0;i3 < 4;i3++){
- unew[i2] += hc->props[os + i3] * u[i3];
+
+ if(hc->rprops[ip1] >= rbound_kludge){
+ //
+ // PROPAGATE U TO NEXT RADIUS IN RPROPS
+ //
+ for(os=pos1 + i*16,i2=0;i2 < 4;i2++,os += 4){
+ unew[i2] = 0.0;
+ for(i3=0;i3 < 4;i3++){
+ unew[i2] += hc->props[os + i3] * u[i3];
+ }
}
- }
- hc_a_equals_b_vector(u,unew,4);
- //
- // PROPAGATE POTEN TO NEXT RADIUS
- //
- os = pos2 + i * 4;
- potnew[0] = poten[0] * hc->ppots[os+0] +
- poten[1] * hc->ppots[os+1];
- poten[1] = poten[0] * hc->ppots[os+2] +
- poten[1] * hc->ppots[os+3];
- poten[0] = potnew[0];
- if(ibv == 0){
+ hc_a_equals_b_vector(u,unew,4);
//
- // ADD DEN * B, WHERE DEN = 0 FOR NO DENSITY CONTRAST
+ // PROPAGATE POTEN TO NEXT RADIUS
//
- dadd = hc->den[i] * b[ninho];
- u[2] += dadd; /* this would have a factor
+ os = pos2 + i * 4;
+ potnew[0] = poten[0] * hc->ppots[os+0] +
+ poten[1] * hc->ppots[os+1];
+ poten[1] = poten[0] * hc->ppots[os+2] +
+ poten[1] * hc->ppots[os+3];
+ poten[0] = potnew[0];
+ if(ibv == 0){
+ //
+ // ADD DEN * B, WHERE DEN = 0 FOR NO DENSITY CONTRAST
+ //
+ dadd = hc->den[i] * b[ninho];
+ u[2] += dadd; /* this would have a factor
grav(i)/hc->grav
*/
+ //
+ // ADD DEN * BETA * B * RDEN
+ //
+ // fprintf(stderr,"%15.5e %15.5e %15.5e %15.5e\n",
+ // beta, hc->den[i], b[ninho],hc->rden[ninho]);
+ poten[1] += hc->psp.beta * dadd * hc->rden[ninho];
+ }
//
- // ADD DEN * BETA * B * RDEN
+ // Changes due to radial density variations
+ //
+ drho = hc->rho_zero[i] - hc->rho_zero[ip1];
+
+ du1 = u[0] * drho/hc->rho_zero[ip1];
+ du2 = du1 * (hc->pvisc[i]+hc->pvisc[ip1]);
+ u[0] += du1;
+ u[2] -= 2.0 * du2 + drho * poten[0];
+ u[3] += du2;
+ //hc_print_vector(poten,2,stderr);
+ //hc_print_vector(u,4,stderr);
//
- // fprintf(stderr,"%15.5e %15.5e %15.5e %15.5e\n",
- // beta, hc->den[i], b[ninho],hc->rden[ninho]);
- poten[1] += hc->psp.beta * dadd * hc->rden[ninho];
- }
- //
- // Changes due to radial density variations
- //
- drho = hc->rho_zero[i] - hc->rho_zero[ip1];
-
- du1 = u[0] * drho/hc->rho_zero[ip1];
- du2 = du1 * (hc->pvisc[i]+hc->pvisc[ip1]);
- u[0] += du1;
- u[2] -= 2.0 * du2 + drho * poten[0];
- u[3] += du2;
- //hc_print_vector(poten,2,stderr);
- //hc_print_vector(u,4,stderr);
- //
- // effects of phase boundary deflections
- //
- if ((jpb < npb)&&(hc->rprops[ip1] > rpb[jpb] - 0.0001)){
- if (ibv == 0) {
- u[2] -= fpb[jpb] * b[ninho];
- poten[1] -= hc->psp.beta * rpb[jpb] * fpb[jpb] * b[ninho] * hc->psp.rho_scale;
+ // effects of phase boundary deflections
+ //
+ if ((jpb < npb)&&(hc->rprops[ip1] > rpb[jpb] - 0.0001)){
+ if (ibv == 0) {
+ u[2] -= fpb[jpb] * b[ninho];
+ poten[1] -= hc->psp.beta * rpb[jpb] * fpb[jpb] * b[ninho] * hc->psp.rho_scale;
+ }
+ jpb++;
}
- jpb++;
+ /* end of l-dependent solver kludge branch */
+ }else{
+ if((verbose)&&(!kludge_warned)){
+ fprintf(stderr,"hc_polsol: applying CMB fixed kludge above %i and shifting CMB to %g km depth for l %i\n",
+ hc->psp.solver_kludge_l,
+ (double)HC_Z_DEPTH(rbound_kludge),l);
+ kludge_warned = TRUE;
+ }
}
+
//
// IF AT A DENSITY CONTRAST, INCREMENT NINHO FOR NEXT ONE
- //
+
if(fabs(hc->den[i]) > HC_EPS_PREC)
ninho++;
//
Modified: mc/1D/hc/trunk/make_tar
===================================================================
--- mc/1D/hc/trunk/make_tar 2012-07-16 12:55:19 UTC (rev 20527)
+++ mc/1D/hc/trunk/make_tar 2012-07-18 22:19:20 UTC (rev 20528)
@@ -2,7 +2,7 @@
#
# create a tar file
#
-ver=${1-1.0.1}
+ver=${1-1.0.2}
date=`date '+%m%d%y'`
tf=$HOME/tmp/hc-$ver.$date.tgz
More information about the CIG-COMMITS
mailing list