From 21ca027a699e61123a927620f29504a3bb415796 Mon Sep 17 00:00:00 2001 From: knarf Date: Tue, 14 Sep 2010 19:45:07 +0000 Subject: Simplify and correct atmosphere handling - This will be taken out of this thorn entirely later, this is only a cleanup to prepare this - This actually also fixes a bug with the atmosphere of multiple stars using the average method, thus the change to the testsuite Remove some unused function arguments of internal function git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TOVSolver/trunk@119 1bdb13ef-5d69-4035-bb54-08abeb3aa7f1 --- src/external.inc | 9 ++--- src/tov.c | 48 ++++++++++--------------- test/test_two_av/dens_3D_diagonal.xg | 66 +++++++++++++++++------------------ test/test_two_av/dens_minimum.xg | 2 +- test/test_two_av/dens_norm1.xg | 2 +- test/test_two_av/dens_norm2.xg | 2 +- test/test_two_av/dens_x_[16][16].xg | 44 +++++++++++------------ test/test_two_av/dens_y_[16][16].xg | 66 +++++++++++++++++------------------ test/test_two_av/dens_z_[16][16].xg | 66 +++++++++++++++++------------------ test/test_two_av/eps_3D_diagonal.xg | 66 +++++++++++++++++------------------ test/test_two_av/eps_minimum.xg | 2 +- test/test_two_av/eps_norm1.xg | 2 +- test/test_two_av/eps_norm2.xg | 2 +- test/test_two_av/eps_x_[16][16].xg | 44 +++++++++++------------ test/test_two_av/eps_y_[16][16].xg | 66 +++++++++++++++++------------------ test/test_two_av/eps_z_[16][16].xg | 66 +++++++++++++++++------------------ test/test_two_av/ham_3D_diagonal.xg | 62 ++++++++++++++++---------------- test/test_two_av/ham_maximum.xg | 2 +- test/test_two_av/ham_norm1.xg | 2 +- test/test_two_av/ham_norm2.xg | 2 +- test/test_two_av/ham_x_[16][16].xg | 52 +++++++++++++-------------- test/test_two_av/ham_y_[16][16].xg | 62 ++++++++++++++++---------------- test/test_two_av/ham_z_[16][16].xg | 62 ++++++++++++++++---------------- test/test_two_av/press_3D_diagonal.xg | 66 +++++++++++++++++------------------ test/test_two_av/press_minimum.xg | 2 +- test/test_two_av/press_norm1.xg | 2 +- test/test_two_av/press_norm2.xg | 2 +- test/test_two_av/press_x_[16][16].xg | 44 +++++++++++------------ test/test_two_av/press_y_[16][16].xg | 66 +++++++++++++++++------------------ test/test_two_av/press_z_[16][16].xg | 66 +++++++++++++++++------------------ test/test_two_av/rho_3D_diagonal.xg | 66 +++++++++++++++++------------------ test/test_two_av/rho_minimum.xg | 2 +- test/test_two_av/rho_norm1.xg | 2 +- test/test_two_av/rho_norm2.xg | 2 +- test/test_two_av/rho_x_[16][16].xg | 44 +++++++++++------------ test/test_two_av/rho_y_[16][16].xg | 66 +++++++++++++++++------------------ test/test_two_av/rho_z_[16][16].xg | 66 +++++++++++++++++------------------ test/test_two_av/tau_3D_diagonal.xg | 66 +++++++++++++++++------------------ test/test_two_av/tau_minimum.xg | 2 +- test/test_two_av/tau_norm1.xg | 2 +- test/test_two_av/tau_norm2.xg | 2 +- test/test_two_av/tau_x_[16][16].xg | 44 +++++++++++------------ test/test_two_av/tau_y_[16][16].xg | 66 +++++++++++++++++------------------ test/test_two_av/tau_z_[16][16].xg | 66 +++++++++++++++++------------------ 44 files changed, 764 insertions(+), 777 deletions(-) diff --git a/src/external.inc b/src/external.inc index 0ab7550..6461c65 100644 --- a/src/external.inc +++ b/src/external.inc @@ -85,8 +85,7 @@ CCTK_INT TOV_Set_Rho_ADM(CCTK_POINTER_TO_CONST cctkGH, &(TOV_press_1d[star_i]), &(TOV_phi_1d[star_i]), &(TOV_rbar_1d[star_i]), &(TOV_r_1d[star_i]), &r_to_star, TOV_Surface[star], - &press, &phi, &r, - 0, 1, NULL); + &press, &phi, &r); press = (press > 0.0) ? press : 0.0; rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]); eps = (rho > 0.0) ? press/(TOV_Gamma[star] - 1.0) / rho : 0.0; @@ -156,8 +155,7 @@ CCTK_INT TOV_Set_Momentum_Source( &(TOV_press_1d[star_i]), &(TOV_phi_1d[star_i]), &(TOV_rbar_1d[star_i]), &(TOV_r_1d[star_i]), &r_to_star, TOV_Surface[star], - &press, &phi, &r, - 0, 1, NULL); + &press, &phi, &r); rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]); gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30); psip = pow(gxx, TOV_Momentum_Psi_Power/4.); @@ -229,8 +227,7 @@ CCTK_INT TOV_Set_Initial_Guess_for_u( &(TOV_press_1d[star_i]), &(TOV_phi_1d[star_i]), &(TOV_rbar_1d[star_i]), &(TOV_r_1d[star_i]), &r_to_star, TOV_Surface[star], - &press, &phi, &r, - 0, 1, NULL); + &press, &phi, &r); rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]); if (rho>0.0) eps = press/(TOV_Gamma[star] - 1.0) / rho; diff --git a/src/tov.c b/src/tov.c index 3c82c44..68e3171 100644 --- a/src/tov.c +++ b/src/tov.c @@ -332,7 +332,6 @@ void TOV_C_Integrate_RHS(CCTK_ARGUMENTS) local_rho = pow(TOV_press_1d[i+1] / TOV_K[star], 1.0 / TOV_Gamma[star]); /* scan for the surface */ - /* if ( (local_rho <= TOV_Atmosphere[star]) ||*/ if ( (local_rho <= 0.0) || (TOV_press_1d[i+1] <= 0.0) ) { @@ -361,8 +360,7 @@ void TOV_C_Integrate_RHS(CCTK_ARGUMENTS) /* match to Schwarzschield */ if (i > TOV_Surface_Index) { - TOV_press_1d[i] = TOV_K[star] * pow(TOV_Atmosphere[star], - TOV_Gamma[star]); + TOV_press_1d[i] = 0.0; TOV_rbar_1d [i] = 0.5 * (sqrt(TOV_r_1d[i]*(TOV_r_1d[i] - 2.0*Surface_Mass)) + TOV_r_1d[i] - Surface_Mass); @@ -444,10 +442,7 @@ void TOV_C_interp_tov_isotropic( CCTK_REAL surface, CCTK_REAL *press_point, CCTK_REAL *phi_point, - CCTK_REAL *r_point, - CCTK_INT TOV_warn_small_grid, - CCTK_INT TOV_zero_atmosphere, - CCTK_REAL *TOV_Atmosphere) + CCTK_REAL *r_point) { DECLARE_CCTK_PARAMETERS CCTK_INT left_index; @@ -459,17 +454,8 @@ void TOV_C_interp_tov_isotropic( *r=TOV_rbar_1d_local[1]; if (*r > TOV_rbar_1d_local[TOV_Num_Radial-2]) { - if (TOV_warn_small_grid) - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Grid not large enough - last point: %.20f.", - TOV_rbar_1d_local[TOV_Num_Radial-1]); - else { - if (!TOV_zero_atmosphere) - *press_point= TOV_K[star] * pow(TOV_Atmosphere[star], - TOV_Gamma[star]); - else - *press_point=0.0; + *press_point= 0.0; M = 0.5 * TOV_r_1d_local[TOV_Num_Radial-1] * (1.0 - exp(2.0*TOV_phi_1d_local[TOV_Num_Radial-1])); *r_point=(2* *r+M)*(2* *r+M)*0.25/ *r; @@ -495,8 +481,7 @@ void TOV_C_interp_tov_isotropic( h * TOV_r_1d_local[left_index+1]; *phi_point = (1.0 - h) * TOV_phi_1d_local[left_index] + h * TOV_phi_1d_local[left_index+1]; - if (!TOV_zero_atmosphere || - (*r_point < surface)) + if (*r_point < surface) *press_point = (1.0 - h) * TOV_press_1d_local[left_index] + h * TOV_press_1d_local[left_index+1]; else @@ -525,7 +510,7 @@ void TOV_C_Exact(CCTK_ARGUMENTS) CCTK_INT LSH_MAX_I; - CCTK_INT i,j,k, i3D, star, star_i; + CCTK_INT i,j,k, i3D, star; CCTK_REAL *r_to_star; CCTK_REAL g_diag, max_g_diag, max_rho, det, sqrt_det; CCTK_REAL my_velx, my_vely, my_velz, my_psi4; @@ -646,7 +631,7 @@ void TOV_C_Exact(CCTK_ARGUMENTS) (y[i3D]-TOV_Position_y[star]) + (z[i3D]-TOV_Position_z[star]) * (z[i3D]-TOV_Position_z[star]) ); - star_i = star * TOV_Num_Radial; + int star_i = star * TOV_Num_Radial; /* do the actual interpolation */ TOV_C_interp_tov_isotropic(star, @@ -654,8 +639,7 @@ void TOV_C_Exact(CCTK_ARGUMENTS) &(TOV_rbar_1d[star_i]), &(TOV_r_1d[star_i]), &(r_to_star[star]), TOV_Surface[star], &(press_point[star]), - &(phi_point[star]), &(r_point[star]), - 0, 0, TOV_Atmosphere); + &(phi_point[star]), &(r_point[star])); /* is some perturbation wanted? */ if (Perturb[star] == 0) @@ -668,8 +652,11 @@ void TOV_C_Exact(CCTK_ARGUMENTS) Pert_Amplitude[star] * cos(PI/2.0 * r[i3D] / TOV_R_Surface[star])); - eps_point[star] = press_point[star] / (TOV_Gamma[star] - 1.0) - / rho_point[star]; + if (rho_point[star] > local_tiny) + eps_point[star] = press_point[star] / (TOV_Gamma[star] - 1.0) + / rho_point[star]; + else + eps_point[star] = 0.0; mu_point[star] = rho_point[star] * (1.0 + eps_point[star]); } /* find out from which star we want to have the data */ @@ -679,7 +666,7 @@ void TOV_C_Exact(CCTK_ARGUMENTS) star=0; max_g_diag = 0.0; max_rho = rho_point[0]; - for (star_i=0; star_i