From 1d77b64f4ef1bd2ad5a28ee48768e6ea12e23e8b Mon Sep 17 00:00:00 2001 From: knarf Date: Mon, 1 Mar 2010 17:52:34 +0000 Subject: Make my old interface from the TOVSolver to TwoPunctures public git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TOVSolver/trunk@79 1bdb13ef-5d69-4035-bb54-08abeb3aa7f1 --- interface.ccl | 56 +++++ param.ccl | 4 + schedule.ccl | 9 + src/external.inc | 642 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 711 insertions(+) diff --git a/interface.ccl b/interface.ccl index ce1f981..ccd7415 100644 --- a/interface.ccl +++ b/interface.ccl @@ -12,3 +12,59 @@ void FUNCTION SpatialDet(CCTK_REAL IN gxx, CCTK_REAL IN gxy, \ CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ CCTK_REAL OUT det) USES FUNCTION SpatialDet + +CCTK_INT FUNCTION Set_Rho_ADM( \ + CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN size, \ + CCTK_REAL ARRAY OUT source, \ + CCTK_REAL ARRAY IN x, \ + CCTK_REAL ARRAY IN y, \ + CCTK_REAL ARRAY IN z \ + ) +PROVIDES FUNCTION Set_Rho_ADM \ + WITH TOV_Set_Rho_ADM \ + LANGUAGE C + +CCTK_INT FUNCTION Set_Momentum_Source( \ + CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN dir, \ + CCTK_INT IN size, \ + CCTK_REAL ARRAY OUT source, \ + CCTK_REAL ARRAY IN x, \ + CCTK_REAL ARRAY IN y, \ + CCTK_REAL ARRAY IN z \ + ) +PROVIDES FUNCTION Set_Momentum_Source \ + WITH TOV_Set_Momentum_Source \ + LANGUAGE C + +CCTK_INT FUNCTION Set_Initial_Guess_for_u( \ + CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN size, \ + CCTK_REAL ARRAY OUT u, \ + CCTK_REAL ARRAY IN x, \ + CCTK_REAL ARRAY IN y, \ + CCTK_REAL ARRAY IN z \ + ) +PROVIDES FUNCTION Set_Initial_Guess_for_u \ + WITH TOV_Set_Initial_Guess_for_u \ + LANGUAGE C + +CCTK_INT FUNCTION Rescale_Sources( \ + CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN size, \ + CCTK_REAL ARRAY IN x, \ + CCTK_REAL ARRAY IN y, \ + CCTK_REAL ARRAY IN z, \ + CCTK_REAL ARRAY IN psi, \ + CCTK_REAL ARRAY IN gxx, \ + CCTK_REAL ARRAY IN gyy, \ + CCTK_REAL ARRAY IN gzz, \ + CCTK_REAL ARRAY IN gxy, \ + CCTK_REAL ARRAY IN gxz, \ + CCTK_REAL ARRAY IN gyz \ + ) +PROVIDES FUNCTION Rescale_Sources \ + WITH TOV_Rescale_Sources \ + LANGUAGE C + diff --git a/param.ccl b/param.ccl index 55602b3..77dc77f 100644 --- a/param.ccl +++ b/param.ccl @@ -119,6 +119,10 @@ CCTK_INT TOV_fake_evolution "Fake evolution by setting ID at every step" STEERAB *:* :: "anything, 0 as off (default), everything else as on" } 0 +STRING TOV_save_to_datafile "Only save data to file and exit" +{ + ".*" :: "Any filename, not used if empty" +} "" shares:HydroBase diff --git a/schedule.ccl b/schedule.ccl index f4c2446..3de374f 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -42,6 +42,15 @@ schedule TOV_C_Integrate_RHS IN TOV_Initial_Data OPTIONS: GLOBAL } "Integrate the 1d equations for the TOV star" +if (!CCTK_Equals(TOV_save_to_datafile,"")) +{ + schedule TOV_write_1D_datafile IN TOV_Initial_Data AFTER TOV_C_Integrate_RHS BEFORE TOV_C_Exact + { + LANG: C + OPTIONS: GLOBAL + } "Save data to file and exit" +} + if (CCTK_Equals(initial_data, "tov") || CCTK_Equals(initial_hydro, "tov") || (TOV_Use_Old_Initial_Data > 0) || TOV_Enforce_Interpolation) { diff --git a/src/external.inc b/src/external.inc index 156c758..cc49146 100644 --- a/src/external.inc +++ b/src/external.inc @@ -1,6 +1,648 @@ +/* for debug output the location of the points */ +void TOV_debug_input_points(CCTK_INT size, + const CCTK_REAL *x, + const CCTK_REAL *y, + const CCTK_REAL *z) +{ + FILE *f; + int i; + f=fopen("test.out", "w"); + for (i=0; i 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; + gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30); + w_lorentz_2 = 1. - gxx*TOV_Velocity_x[star]*TOV_Velocity_x[star] + - gxx*TOV_Velocity_y[star]*TOV_Velocity_y[star] + - gxx*TOV_Velocity_z[star]*TOV_Velocity_z[star]; + if (rho > 0.0) + /* source[i]=gxx*gxx*rho*(1.0+eps); */ + source[i]=gxx*gxx* (w_lorentz_2*rho + + (w_lorentz_2*TOV_Gamma[star]/(TOV_Gamma[star]-1.)-1.)* + press); + else + source[i]=0.0; + } + return 1; +} + +/* This function provides an interface for other thorns which can input a list + * of points (p_x, p_y, p_z) with a size of 'size' and get back the momentum + * sources at these points (source) */ +/* NOTE: when used with two or more NS, there might be problems if they are + * very close to each other, because this function does not use the fancy + * algorithms above to determine from which NS the data for a perticular grid + * point should be taken, but just uses the closest one */ +CCTK_INT TOV_Set_Momentum_Source( + CCTK_POINTER_TO_CONST cctkGH, + CCTK_INT dir, + CCTK_INT size, + CCTK_REAL *source, + const CCTK_REAL *x, + const CCTK_REAL *y, + const CCTK_REAL *z ) +{ + DECLARE_CCTK_PARAMETERS + CCTK_REAL r_to_star, press, phi, r, gxx, rho, rho_eps, psip, w_lorentz_2; + CCTK_INT star_i=0, star=0; + CCTK_INT i; + + assert(TOV_Surface!=0); + assert(TOV_R_Surface!=0); + assert(TOV_Atmosphere!=0); + + assert(TOV_r_1d!=0); + assert(TOV_rbar_1d!=0); + assert(TOV_press_1d!=0); + assert(TOV_phi_1d!=0); + assert(TOV_m_1d!=0); + + /* CCTK_INFO("Request for Momentum-Sources from TOVSolverC"); */ + /* loop over all points we got */ + + for (i=0; i0.0) + eps = press/(TOV_Gamma[star] - 1.0) / rho; + else + eps = 0.0; + gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30); + u[i]=pow(gxx, 0.25) - 1.0; + } + return 1; +} + + +void bisection (CCTK_REAL *vx, CCTK_REAL *vy, CCTK_REAL *vz, + CCTK_REAL *rhoold, CCTK_REAL *rhonew, + CCTK_REAL *w_lorentz_2, CCTK_REAL *v_2, + CCTK_REAL Gamma, CCTK_REAL K, + CCTK_REAL source, + CCTK_REAL mom_source_x, + CCTK_REAL mom_source_y, + CCTK_REAL mom_source_z, + CCTK_REAL my_psi4, + CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL rho_orig) +{ + CCTK_REAL lb, ub, f; + lb=0.0; + ub=2*rho_orig; + do + { + *rhonew = (ub-lb)/2.; + *vx = mom_source_x/my_psi4/ + (source/my_psi4/my_psi4+K * pow(*rhonew, Gamma)); + *vy = mom_source_y/my_psi4/ + (source/my_psi4/my_psi4+K * pow(*rhonew, Gamma)); + *vz = mom_source_z/my_psi4/ + (source/my_psi4/my_psi4+K * pow(*rhonew, Gamma)); + *v_2 = my_psi4*(*vx**vx+*vy**vy+*vz**vz); + *w_lorentz_2 = 1.- *v_2; + f = my_psi4*my_psi4*( + *w_lorentz_2 * *rhonew + + (*w_lorentz_2*Gamma/(Gamma-1.)-1.) * + K * pow(*rhonew, Gamma)) - source; + if (f < 0.0) + ub = *rhonew; + else + lb = *rhonew; + *rhoold = *rhonew; + } while (ub-lb < 1.e-10); +} + +void newton_raphson(CCTK_REAL *vx, CCTK_REAL *vy, CCTK_REAL *vz, + CCTK_REAL *rhoold, CCTK_REAL *rhonew, + CCTK_REAL *w_lorentz_2, CCTK_REAL *v_2, + CCTK_REAL Gamma, CCTK_REAL K, + CCTK_REAL source, + CCTK_REAL mom_source_x, + CCTK_REAL mom_source_y, + CCTK_REAL mom_source_z, + CCTK_REAL my_psi4, + CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL rho_orig) +{ + int count=0; + CCTK_REAL d_w_lorentz_2, f, df; + do + { + count++; + *vx = mom_source_x/my_psi4/ + (source/my_psi4/my_psi4+K * pow(*rhonew, Gamma)); + *vy = mom_source_y/my_psi4/ + (source/my_psi4/my_psi4+K * pow(*rhonew, Gamma)); + *vz = mom_source_z/my_psi4/ + (source/my_psi4/my_psi4+K * pow(*rhonew, Gamma)); + *v_2 = my_psi4*(*vx**vx+*vy**vy+*vz**vz); + if (count > 100) + { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, \ + "Failed to converge IVP rescaling\nat (%.3g %.3g %.3g)\n"\ + " orig_rho:%g rho:%g, v^2=%g, rel_diff=%g", \ + x, y, z, rho_orig, *rhonew, *v_2, + fabs(*rhonew-*rhoold)/fabs(*rhonew)); + CCTK_WARN(1, "trying slow bisection method"); + bisection(vx,vy,vz,rhoold,rhonew,w_lorentz_2,v_2, + Gamma, K, + source, mom_source_x, mom_source_y, mom_source_z, + my_psi4, x, y, z, rho_orig); + return; + } + *w_lorentz_2 = 1.- *v_2; + f = my_psi4*my_psi4*( + *w_lorentz_2 * *rhonew + + (*w_lorentz_2*Gamma/(Gamma-1.)-1.) * + K * pow(*rhonew, Gamma)) - source; + d_w_lorentz_2 = + 2. * (*w_lorentz_2 - 1.) / + (source/my_psi4/my_psi4+K*pow(*rhonew, Gamma)) * + K * Gamma * pow(*rhonew, Gamma-1.); + df = my_psi4*my_psi4*( + d_w_lorentz_2**rhonew + *w_lorentz_2 + + d_w_lorentz_2*Gamma/(Gamma-1.)* + K * pow(*rhonew, Gamma) + + (*w_lorentz_2*Gamma/(Gamma-1.)-1.) * + K * Gamma * pow(*rhonew, Gamma-1.)); + *rhoold = *rhonew; + *rhonew = *rhoold - f/df; + } + while ( (fabs(*rhonew - *rhoold)/fabs(*rhonew) > 1.e-12) && + (fabs(*rhonew - *rhoold) > 1.e-12) + ); +} + +CCTK_INT TOV_Rescale_Sources( + CCTK_POINTER_TO_CONST cctkGH, + CCTK_INT size, + const CCTK_REAL *x, + const CCTK_REAL *y, + const CCTK_REAL *z, + const CCTK_REAL *psi, + const CCTK_REAL *gxx, + const CCTK_REAL *gyy, + const CCTK_REAL *gzz, + const CCTK_REAL *gxy, + const CCTK_REAL *gxz, + const CCTK_REAL *gyz) +{ + DECLARE_CCTK_PARAMETERS + CCTK_REAL *rho, *press, *eps, *dens, *tau, *w_lorentz, *vel0, *vel1, *vel2, + *scon0, *scon1, *scon2; + CCTK_REAL *rho_p, *press_p, *eps_p, *dens_p, *tau_p, *w_lorentz_p, *vel0_p, + *vel1_p, *vel2_p, *scon0_p, *scon1_p, *scon2_p; + CCTK_REAL *rho_p_p, *press_p_p, *eps_p_p, *dens_p_p, *tau_p_p, *w_lorentz_p_p, + *vel0_p_p, *vel1_p_p, *vel2_p_p, *scon0_p_p, *scon1_p_p, *scon2_p_p; + CCTK_REAL *u, *source, *mom_source[3]; + CCTK_REAL psip, my_psi4, sqrt_det, vx, vy, vz, v_2, + w_lorentz_2, d_w_lorentz_2, D_h_w; + CCTK_REAL vlowx, vlowy, vlowz; + CCTK_INT i, type; + CCTK_REAL rhoold, rhonew, f, df; + int i3D; + static int debug = 1; + FILE *debugfile = NULL; + + CCTK_INFO("Rescaling Sources"); + + /* get Whisky variables */ + rho =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "HydroBase::rho"); + press =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "HydroBase::press"); + eps =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "HydroBase::eps"); + dens =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "Whisky::dens"); + tau =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "Whisky::tau"); + w_lorentz=(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "Whisky::w_lorentz"); + vel0 =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "HydroBase::vel[0]"); + vel1 =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "HydroBase::vel[1]"); + vel2 =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "HydroBase::vel[2]"); + scon0 =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "Whisky::scon[0]"); + scon1 =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "Whisky::scon[1]"); + scon2 =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "Whisky::scon[2]"); + + rho_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "HydroBase::rho"); + press_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "HydroBase::press"); + eps_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "HydroBase::eps"); + dens_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "Whisky::dens"); + tau_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "Whisky::tau"); + w_lorentz_p=(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "Whisky::w_lorentz"); + vel0_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "HydroBase::vel[0]"); + vel1_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "HydroBase::vel[1]"); + vel2_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "HydroBase::vel[2]"); + scon0_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "Whisky::scon[0]"); + scon1_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "Whisky::scon[1]"); + scon2_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "Whisky::scon[2]"); + + rho_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "HydroBase::rho"); + press_p_p=(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "HydroBase::press"); + eps_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "HydroBase::eps"); + dens_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "Whisky::dens"); + tau_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "Whisky::tau"); + w_lorentz_p_p=(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "Whisky::w_lorentz"); + vel0_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "HydroBase::vel[0]"); + vel1_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "HydroBase::vel[1]"); + vel2_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "HydroBase::vel[2]"); + scon0_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "Whisky::scon[0]"); + scon1_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "Whisky::scon[1]"); + scon2_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "Whisky::scon[2]"); + + /* get u */ + u=calloc(size, sizeof(CCTK_REAL)); + TOV_Set_Initial_Guess_for_u(cctkGH, size, u, x, y, z); + + /* get the old source */ + source=calloc(size, sizeof(CCTK_REAL)); + TOV_Set_Rho_ADM(cctkGH, size, source, x, y, z); + mom_source[0]=calloc(size, sizeof(CCTK_REAL)); + mom_source[1]=calloc(size, sizeof(CCTK_REAL)); + mom_source[2]=calloc(size, sizeof(CCTK_REAL)); + TOV_Set_Momentum_Source(cctkGH, 0, size, mom_source[0], x, y, z); + TOV_Set_Momentum_Source(cctkGH, 1, size, mom_source[1], x, y, z); + TOV_Set_Momentum_Source(cctkGH, 2, size, mom_source[2], x, y, z); + + if (debug) + { + debugfile=fopen("tovdebug.dat", "w"); + } + /* loop over all points we got */ + for (i3D=0; i3D0.0) + { + int count=0; + if (debug && (fabs(z[i3D])<1.e-15)) + { + fprintf(debugfile, + "%.15g %.15g %.15g %.15g %.15g %.15g ", + x[i3D], y[i3D], source[i3D], + mom_source[0][i3D], mom_source[1][i3D], mom_source[1][i3D]); + } + if (psi != NULL) + my_psi4 = pow(psi[i3D], 4.0)*gxx[i3D]; + else + my_psi4 = gxx[i3D]; + /* we do not need u, we need gxx here, which will be stored in u! */ + u[i3D] = pow(u[i3D]+1.0, 4.0); + /* We also store the original J^i in mom_source[i] */ + psip = pow(u[i3D], TOV_Momentum_Psi_Power/4.); + mom_source[0][i3D] /= psip; + mom_source[1][i3D] /= psip; + mom_source[2][i3D] /= psip; + + if (rho[i3D] < 1.e-10) + rho[i3D] = 1.e-10; + rhoold = rhonew = rho[i3D]; +#ifdef OLD_NR + do + { + count++; + vx = mom_source[0][i3D]/my_psi4/ + (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0])); + vy = mom_source[1][i3D]/my_psi4/ + (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0])); + vz = mom_source[2][i3D]/my_psi4/ + (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0])); + v_2 = my_psi4*(vx*vx+vy*vy+vz*vz); + if (count > 100) + CCTK_VWarn(count<110, __LINE__, __FILE__, CCTK_THORNSTRING, \ + "Failed to converge IVP rescaling\nat (%.3g %.3g %.3g)\n"\ + " orig_rho:%g rho:%g, v^2=%g, rel_diff=%g", \ + x[i3D], y[i3D], z[i3D], rho[i3D], rhonew, v_2, + fabs(rhonew-rhoold)/fabs(rhonew)); + w_lorentz_2 = 1.- v_2; + /* + f = TOV_K[0] / (TOV_Gamma[0]-1.0) * pow(rhonew, TOV_Gamma[0]) + + rhonew - source[i3D]; + df = TOV_K[0] * TOV_Gamma[0] / (TOV_Gamma[0]-1.0) * + pow(rhonew, TOV_Gamma[0]-1.0) + 1.0; + */ + f = my_psi4*my_psi4*( + w_lorentz_2 * rhonew + + (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) * + TOV_K[0] * pow(rhonew, TOV_Gamma[0])) - source[i3D]; + /* + df= w_lorentz_2 + + (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) * + TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.); + */ + d_w_lorentz_2 = + 2. * (w_lorentz_2 - 1.) / + (source[i3D]/my_psi4/my_psi4+TOV_K[0]*pow(rhonew, TOV_Gamma[0])) * + TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.); + df = my_psi4*my_psi4*( + d_w_lorentz_2*rhonew + w_lorentz_2 + + d_w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)* + TOV_K[0] * pow(rhonew, TOV_Gamma[0]) + + (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) * + TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.)); + /* + df= 1. + my_psi4*v_2/(rhonew*rhonew) + + ((TOV_Gamma[0]+(2.-TOV_Gamma[0])*my_psi4*v_2/(rhonew*rhonew)) + / (TOV_Gamma[0]-1.) - 1.) + * TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.); + */ + + rhoold = rhonew; + rhonew = rhoold - f/df; + } + while ( (fabs(rhonew-rhoold)/fabs(rhonew) > 1.e-12) && + (fabs(rhonew - rhoold) > 1.e-12) + ); + if (count==1) + printf("Fast rescale! %g %g %g %g\n", u[i3D], my_psi4, + mom_source[0][i3D]/u[i3D]/ + (source[i3D]/u[i3D]/u[i3D]+TOV_K[0]* + pow(rhonew, TOV_Gamma[0])), + mom_source[0][i3D]/my_psi4/ + (source[i3D]/u[i3D]/my_psi4+TOV_K[0]* + pow(rhonew, TOV_Gamma[0]))); +#else + newton_raphson(&vx,&vy,&vz,&rhoold,&rhonew,&w_lorentz_2,&v_2, + TOV_Gamma[0], TOV_K[0], + source[i3D], + mom_source[0][i3D],mom_source[1][i3D],mom_source[2][i3D], + my_psi4, x[i3D], y[i3D], z[i3D], rho[i3D]); +#endif + rho[i3D] = rhonew; + press[i3D] = TOV_K[0] * pow(rhonew, TOV_Gamma[0]); + eps[i3D] = TOV_K[0] * pow(rhonew, TOV_Gamma[0] - 1.0) / + (TOV_Gamma[0] - 1.0); + sqrt_det = pow(my_psi4, 1.5); + w_lorentz[i3D] = sqrt(w_lorentz_2); + dens[i3D] = sqrt_det * w_lorentz[i3D] * rho[i3D]; + tau[i3D] = sqrt_det * (rho[i3D] *w_lorentz[i3D]*w_lorentz[i3D] * + (1.0 + eps[i3D]) + + press[i3D]*(w_lorentz[i3D]*w_lorentz[i3D]-1.0) ) - + dens[i3D]; + vel0[i3D] = vx; + vel1[i3D] = vy; + vel2[i3D] = vz; + D_h_w = dens[i3D] * (1 + eps[i3D] + press[i3D]/rho[i3D]) * w_lorentz[i3D]; + + vlowx = gxx[i3D]*vel0[i3D] + gxy[i3D]*vel1[i3D] + gxz[i3D]*vel2[i3D]; + vlowy = gxy[i3D]*vel0[i3D] + gyy[i3D]*vel1[i3D] + gyz[i3D]*vel2[i3D]; + vlowz = gxz[i3D]*vel0[i3D] + gyz[i3D]*vel1[i3D] + gzz[i3D]*vel2[i3D]; + + scon0[i3D] = D_h_w * vlowx; + scon1[i3D] = D_h_w * vlowy; + scon2[i3D] = D_h_w * vlowz; + + if (debug && (fabs(z[i3D])<1.e-15)) + { + fprintf(debugfile, + "%.15g %.15g %.15g %.15g\n", /*7*/ + gxx[i3D]*gxx[i3D]*( /* ham - source*/ + w_lorentz_2*rho[i3D]+ + (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.)*press[i3D]), + gxx[i3D]*gxx[i3D]*( /* ham - source */ + w_lorentz_2*(rho[i3D]*(1.+eps[i3D])+press[i3D])-press[i3D]), + (w_lorentz_2 * rho[i3D] + /* momx - source */ + w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.) * press[i3D])* + gxx[i3D]*vel0[i3D], + w_lorentz_2 * (rho[i3D]* /* momx - source */ + (1.+eps[i3D])+ press[i3D])* + gxx[i3D]*vel0[i3D] + ); + } + + /* FIXME: check for Atmosphere */ + } + } + if (debug) + fclose(debugfile); + + switch(*(const CCTK_INT*)CCTK_ParameterGet("TOV_Populate_Timelevels", + "Whisky_TOVSolverC", &type)) + { + case 3: + TOV_Copy(size, rho_p_p, rho); + TOV_Copy(size, eps_p_p, eps); + TOV_Copy(size, vel0_p_p, vel0); + TOV_Copy(size, vel1_p_p, vel1); + TOV_Copy(size, vel2_p_p, vel2); + TOV_Copy(size, dens_p_p, dens); + TOV_Copy(size, tau_p_p, tau); + TOV_Copy(size, scon0_p_p, scon0); + TOV_Copy(size, scon1_p_p, scon1); + TOV_Copy(size, scon2_p_p, scon2); + TOV_Copy(size, w_lorentz_p_p, w_lorentz); + case 2: + TOV_Copy(size, rho_p, rho); + TOV_Copy(size, eps_p, eps); + TOV_Copy(size, vel0_p, vel0); + TOV_Copy(size, vel1_p, vel1); + TOV_Copy(size, vel2_p, vel2); + TOV_Copy(size, dens_p, dens); + TOV_Copy(size, tau_p, tau); + TOV_Copy(size, scon0_p, scon0); + TOV_Copy(size, scon1_p, scon1); + TOV_Copy(size, scon2_p, scon2); + TOV_Copy(size, w_lorentz_p, w_lorentz); + } + free(source); + free(u); + return 1; +} + +void TOV_write_1D_datafile(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + int i; + + if (TOV_Num_TOVs > 1) + CCTK_WARN(0, "Writing a data file for multiple stars is not (yet) " + "supported"); + FILE *file; + file = fopen(TOV_save_to_datafile, "w"); + fprintf(file, "TOVSolver data file\n"); + fprintf(file, "version 1.0\n"); + fprintf(file, "TOV_Num_Radial %d\n", TOV_Num_Radial); + fprintf(file, "\n"); + for (i=0; i < TOV_Num_Radial; i++) + { + fprintf(file, "%g %g %g %g %g\n", + TOV_rbar_1d[i], + TOV_r_1d[i], + TOV_m_1d[i], + TOV_phi_1d[i], + TOV_press_1d[i]); + } + fclose(file); + CCTK_WARN(0, "Simulation stopped as requested after writing data to file"); +} + void TOVSolverC_export_local_variables(CCTK_REAL **exported_TOV_press_1d, CCTK_REAL **exported_TOV_m_1d, CCTK_REAL **exported_TOV_phi_1d, CCTK_REAL **exported_TOV_rbar_1d, CCTK_REAL **exported_TOV_r_1d) { + int i; + *exported_TOV_press_1d = TOV_press_1d; *exported_TOV_m_1d = TOV_m_1d; *exported_TOV_phi_1d = TOV_phi_1d; -- cgit v1.2.3