diff options
Diffstat (limited to 'src/GRHydro_Bondi.c')
| -rw-r--r-- | src/GRHydro_Bondi.c | 46 |
1 files changed, 27 insertions, 19 deletions
diff --git a/src/GRHydro_Bondi.c b/src/GRHydro_Bondi.c index e1b3d4c..9fefdf6 100644 --- a/src/GRHydro_Bondi.c +++ b/src/GRHydro_Bondi.c @@ -37,6 +37,7 @@ #include <math.h> #include <stdio.h> #include <stdlib.h> +#include <assert.h> // set to 1 to tracing output #define LTRACE 1 @@ -129,7 +130,7 @@ static void dxc_dxs_ks_calc(CCTK_REAL *x_cart, CCTK_REAL *x_spher, CCTK_REAL dxc for( i = 0 ; i < NDIM ; i++ ) { for( j = 0 ; j < NDIM ; j++ ) { - dxc_dxs[0][i] = 0. ; + dxc_dxs[j][i] = 0. ; } } @@ -339,7 +340,6 @@ static void iso_cart_to_bl_spher_pos(CCTK_REAL *x_cart, CCTK_REAL *x_spher) ****************************************************************************/ static void setutcon(CCTK_REAL *vcon, CCTK_REAL gcov[][NDIM]) { - int i,j; CCTK_REAL d,b,c; d=gcov[TT][TT]; @@ -367,7 +367,7 @@ static void setutcon(CCTK_REAL *vcon, CCTK_REAL gcov[][NDIM]) ****************************************************************************/ static void bl_gcov_func( CCTK_REAL *x, CCTK_REAL gcov[NDIM][NDIM]) { - int i,j,k ; + int i,j ; CCTK_REAL sth,cth,s2,r2,DD,mu ; CCTK_REAL r,th; @@ -422,7 +422,6 @@ static void bl_gcov_func( CCTK_REAL *x, CCTK_REAL gcov[NDIM][NDIM]) static void set_bondi_parameters( CCTK_REAL M_in, CCTK_REAL Mdot_in, CCTK_REAL rs_in, CCTK_REAL gam ) { - CCTK_REAL my_pi, checkp; CCTK_REAL gtemp; /* Set the solution-determining parameters: */ @@ -499,19 +498,18 @@ static int gnr_bondi( CCTK_REAL x[], int n, CCTK_REAL [][NEWT_DIM_B], CCTK_REAL *, CCTK_REAL *, int) ) { - CCTK_REAL f, df, dx[NEWT_DIM_B], x_old[NEWT_DIM_B], resid[NEWT_DIM_B], + CCTK_REAL f, df, dx[NEWT_DIM_B], resid[NEWT_DIM_B], jac[NEWT_DIM_B][NEWT_DIM_B]; - CCTK_REAL errx, x_orig[NEWT_DIM_B]; - int n_iter, id, jd, i_extra, doing_extra; + CCTK_REAL errx; + int n_iter, id, i_extra, doing_extra; - int keep_iterating, i_increase; + int keep_iterating; // Initialize various parameters and variables: errx = 1. ; df = f = 1.; i_extra = doing_extra = 0; - for( id = 0; id < n ; id++) x_old[id] = x_orig[id] = x[id] ; n_iter = 0; @@ -522,11 +520,7 @@ static int gnr_bondi( CCTK_REAL x[], int n, (*funcd) (x, dx, resid, jac, &f, &df, n); /* returns with new dx, f, df */ - /* Save old values before calculating the new: */ errx = 0.; - for( id = 0; id < n ; id++) { - x_old[id] = x[id] ; - } /* don't use line search : */ for( id = 0; id < n ; id++) { @@ -662,18 +656,20 @@ static int find_bondi_solution( CCTK_REAL r, CCTK_REAL *rho, CCTK_REAL *u, CCTK_ if( *rho < 0. ) { if( r > 0.9*rs && r < 1.1*rs ) { - rho_guess = rhos; - *rho = rho_guess; + *rho = rhos; } else { // rhotmp = (sqrt(Qdot) - 1.) * (gamma_eos - 1.) / ( gamma_eos * K ); // rho_guess = pow( rhotmp , (1./(gamma_eos - 1.)) ); if(r < rs) { ur = pow(r,-0.5) ; } else { ur = 0.5*pow(r,-1.5) ; } - rho_guess = Mdot / (4.*M_PI * r * r * ur); - *rho = rho_guess; + *rho = Mdot / (4.*M_PI * r * r * ur); } - } + } + + // safe guess value for multiple tries + rho_guess = *rho; + // set global variables needed by residual function: r_sol = r ; @@ -783,6 +779,10 @@ static void calc_vel_bondi( CCTK_REAL vtmp, CCTK_REAL x[NDIM], CCTK_REAL x_sphe DLOOP2 { ucon[i] += dxc_dxs[i][j] * ucon_bl[j] ; } break; + default: + assert(0 && "Internal error"); + return; /* NOTREACHED */ + } return; @@ -1007,6 +1007,8 @@ void GRHydro_Bondi(CCTK_ARGUMENTS) while( (j < (nbondi-1)) && (bad_bondi[j]) ) { j++; } if( (j==(nbondi-1)) && bad_bondi[j] ) { CCTK_WARN(CCTK_WARN_ABORT,"No available points with which to interpolate!\n"); + ileft = 0; /* NOTREACED */ + iright = 0; /* NOTREACHED */ } else { ileft = j ; @@ -1016,6 +1018,7 @@ void GRHydro_Bondi(CCTK_ARGUMENTS) while( (j < (nbondi-1)) && (bad_bondi[j]) ) { j++; } if( (j==(nbondi-1)) && bad_bondi[j] ) { CCTK_WARN(CCTK_WARN_ABORT,"Need another point to the right but cannot find one! \n"); + iright = 0; /* NOTREACHED */ } else { iright = j ; @@ -1036,6 +1039,7 @@ void GRHydro_Bondi(CCTK_ARGUMENTS) while( (j > 0) && (bad_bondi[j]) ) { j--; } if( (j==0) && bad_bondi[j] ) { CCTK_WARN(CCTK_WARN_ABORT,"Need another point to the left but cannot find one! \n"); + ileft = 0; /* NOTREACED */ } else { ileft = j; @@ -1059,7 +1063,7 @@ void GRHydro_Bondi(CCTK_ARGUMENTS) #if(LTRACE) for(i=0; i < nbondi; i++) { - CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,"radial-bondisoln %12d %26.16e %26.16e %26.16e %26.16e %12d \n",i,r_bondi[i],rho_bondi[i],u_bondi[i],v_bondi[i], bad_bondi[i]); + CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"radial-bondisoln %12d %26.16e %26.16e %26.16e %26.16e \n",i,r_bondi[i],rho_bondi[i],u_bondi[i],v_bondi[i]); } #endif @@ -1089,6 +1093,10 @@ void GRHydro_Bondi(CCTK_ARGUMENTS) iso_cart_to_bl_spher_pos(xpos, x_spher); break; + default: + assert(0 && "Internal error"); + return; /* NOTREACHED */ + } rspher = x_spher[RR]; |