Class to solve the Tolman-Oppenheimer-Volkov equations. More...
#include <tov_solve.h>
See the Solution of the Tolman-Oppenheimer-Volkov Equations section of the User's Guide for the mathematical background.
This class uses adaptive integration to compute the gravitational mass, the radius, the baryonic mass (if the EOS supplies the baryon density), and the gravitational potential (if requested). The equation of state may be changed at any time, by specifying the appropriate eos_tov object
Basic Usage
After specifying the EOS through tov_solve::set_eos(), one calls either tov_solve::mvsr(), tov_solve::fixed(), tov_solve::max() or tov_solve::fixed_pr() to compute the mass versus radius curve, the profile of a star of a given fixed mass, the profile of the maximum mass star, or the profile of a star with a fixed central pressure. These functions all generate results in the form of a table_units object, which can be obtained from tov_solve::get_results().
Screen output:
verbose=0
- Nothingverbose=1
- Basic informationverbose=2
- For each profile computation, report solution information at every kilometer.verbose=3
- Report profile information at every 20 grid points and output the final interpolation to zero pressure. A keypress is required after each profile.Mass versus radius curve
The neutron star mass versus radius curve is constructed by computing profiles of all neutron stars with a range of central pressures. This range is from prbegin to prend, and the ratio between successive central pressures is specified in princ (making a logarithmic central pressure grid).
Profiles of fixed mass
Profiles for a fixed gravitational mass are computed by solving for the correct central pressure. In order to ensure that the solver does not accidentally select a central pressure beyond the maximum mass neutron star, the profile with the maximum mass is computed beforehand automatically (using max()). The intial guess to the solver is always the value of fixed_pr_guess, which defaults to . Alternatively, the user can specify the central pressure of the maximum mass star so that it does not have to be computed.
In order to handle multiply-branched mass-radius relations, the value of fixed_pr_guess must be specified in order to ensure that the correct profile is generated. Even if fixed_pr_guess is specified, there is no guarantee that the Newton-Raphson will select the correct profile (though it is very likely if the guess for the central pressure is sufficiently accurate).
The fixed() function does not support computing profiles with central pressures with gravitational masses larger than the maximum value. For this, use fixed_pr() .
Profile for the maximum mass star
This is provided by the function max() . Internally, this uses the minimizer specified by set_min() or the default minimizer, def_min, to minimize the function max_fun() . In order to generate a good initial guess, the max() function begins by looping over central pressures from max_begin to max_end and choosing the best guess from that set of configurations.
Profile for a fixed central pressure
This is provided by the function fixed_pr() and is relatively fast because it does not require any solving or minimizing. This function allows central pressures larger than that of the maximum mass star.
Output tables
The functions tov_solve::fixed(), tov_solve::fixed_pr(), and tov_solve::max() produce output tables which represent the profile of the neutron star of the requested mass. The columns are
gm
, the enclosed gravitational mass in r
, the radial coordinate in gp
, the gravitational potential (unitless) when calc_gpot is truerjw
, the value of omega_rat
, the value of bm
, the baryonic mass in pr
, the pressure in user-specified unitsed
, the energy density in user-specified unitsnb
, the baryon density in user-specified units (if eos_tov::baryon_column is true)sg
, the local surface gravity (in rs
, the local redshift (unitless),dmdr
, the derivative of the enclosed gravitational mass in dlogpdr
, the derivative of the natural logarithm of the pressuredgpdr
, the derivative of the gravitational potential in dbmdr
, the derivative of the enclosed baryonic mass (if eos_tov::baryon_column is true).The function tov_solve::mvsr() produces a different kind of output table corresponding to the mass versus radius curve. Some points on the curve may correspond to unstable configurations.
gm
, the total gravitational mass in r
, the radius in gp
, the gravitational potential in the center (unitless) when calc_gpot is truerjw
, the value of omega_rat
, the value of bm
, total the baryonic mass in pr
, the central pressure in user-specified unitsed
, the central energy density in user-specified unitsnb
, the central baryon density in user-specified units (if eos_tov::baryon_column is true)sg
, the surface gravity (in rs
, the redshift at the surface,dmdr
, the derivative of the gravitational massdlogpdr
, the derivative of the natural logarithm of the pressuredgpdr
, the derivative of the gravitational potential in dbmdr
, the derivative of the enclosed baryonic mass (if eos_tov::baryon_column is true).Unit systems
By default, this class operates with energy density and pressure in units of and baryon density in units of
.
The function set_units(std::string,std::string,std::string) allows one to use different unit systems for energy density, pressure, and baryon density. The following list of units of energy density and pressure are hard-coded into the library and always work:
"g/cm^3"
,"erg/cm^3"
,"dyne/cm^2"
,"MeV/fm^3"
,"1/fm^4"
, and"Msun/km^3"
(i.e. The list of hard-coded units for baryon density are:
"1/m^3"
,"1/cm^3"
, and"1/fm^3"
.Other units choices will work if the conversion is either already added to the global unit conversion object (from o2scl_settings.get_convert_units()
) or the global unit conversion object is able to compute them by a system()
call to GNU units
(see documentation in convert_units). Note that the choice of what units the tables are produced in is independent of the unit system specified in the associated eos_tov object, i.e. the input EOS and output EOS units need not be the same.
Alternatively, using set_units(double,double,double) allows one to specify the conversion factors directly without using the global unit conversion object.
Accuracy
The present code, as demonstrated in the tests, gives the correct central pressure and energy density of the analytical solution by Buchdahl to within less than 1 part in .
Rotation
Rotation is considered if tov_solve::ang_vel is set to true
. This adds two more differential equations to solve for each central pressure.
The differential equation for (see the section in the User's Guide called Solution of the Tolman-Oppenheimer-Volkov Equations ) is independent of the relative scale for
and
. First, one rescales
and rewrites everything in terms of
and
. Then, pick a central pressure,
, arbitrary values for
and
, and integrate
Afterwards, shift by a constant to ensure the correct value at
, and multiply
by a constant to ensure that
holds for the new potential
. Then, multiply
by a constant to ensure that
The functions and
are stored in columns called
"omega_rat"
and "rjw"
, respectively. One can compute the baryonic mass by integration or by adding one additional differential equation, bringing the total to six.
The moment of inertia is (see Solution of the Tolman-Oppenheimer-Volkov Equations),
where the last fraction is stored in domega_rat . For an object named ts
of type tov_solve after a call to fixed(), max(), or fixed_pr(), the moment of inertia can be computed with, e.g.
Convergence details
By default, if the TOV solver fails to converge, the error handler is called and an exception is thrown. If o2scl::tov_solve::err_nonconv is false, then o2scl::tov_solve::mvsr(), o2scl::tov_solve::fixed(), and o2scl::tov_solve::max(), return an integer which gives some information about why the solver failed to converge.
If err_nonconv is set to false, then the fixed() function temporarily sets the value of both o2scl::mroot::err_nonconv for the current solver and o2scl::jacobian::err_nonconv for the jacobian object, o2scl::mroot_hybrids::def_jac, of def_solver equal to false.
Other details
The ODE solution is stored in a buffer which can be directly accessed using o2scl::tov_solve::get_rkx(), o2scl::tov_solve::get_rky(), and o2scl::tov_solve::get_rkdydx(). In the case that the ODE steps are small enough that there isn't enough space in the buffers, the ODE solution data is thinned by a factor of two to allow for the remaining ODE steps to be stored. The size of the buffers can be controlled in buffer_size .
If o2scl::tov_solve::reformat_results is true (the default), then the results are placed in a shared pointer to the table_units object which can be accessed using o2scl::tov_solve::get_results(). The maximum size of the output table can be controlled with max_table_size. The output table may be smaller than this, as it cannot be larger than the number of steps stored in the buffer.
gsl_efailed
without calling the error handler in the case that the solver can recover gracefully from, for example, a negative pressure.Definition at line 332 of file tov_solve.h.
Public Types | |
typedef boost::numeric::ublas::vector< double > | ubvector |
typedef boost::numeric::ublas::matrix< double > | ubmatrix |
typedef boost::numeric::ublas::matrix_column< ubmatrix > | ubmatrix_column |
typedef boost::numeric::ublas::matrix_row< ubmatrix > | ubmatrix_row |
Public Member Functions | |
Basic operation | |
void | set_eos (eos_tov &ter) |
Set the EOS to use. | |
void | set_units (double s_efactor=1.0, double s_pfactor=1.0, double s_nbfactor=1.0) |
Set output units for the table. | |
void | set_units (std::string eunits="", std::string punits="", std::string nunits="") |
Set output units for the table. | |
virtual int | mvsr () |
Calculate the mass vs. radius curve. | |
virtual int | fixed (double target_mass, double pmax=1.0e20) |
Calculate the profile of a star with fixed mass. More... | |
virtual int | fixed_pr (double pcent) |
Calculate the profile of a star with a specified central pressure. | |
virtual int | max () |
Calculate the profile of the maximum mass star. | |
virtual void | make_table () |
Construct a table from the results. More... | |
std::shared_ptr< table_units<> > | get_results () |
Return the results data table. | |
void | set_table (std::shared_ptr< table_units<> > t) |
Return the results data table. More... | |
Get internally formatted results (in internal unit system) | |
const ubvector & | get_rkx () const |
Get the vector for the radial grid. | |
const std::vector< ubvector > & | get_rky () const |
Get a list of vectors for the ODEs. | |
const std::vector< ubvector > & | get_rkdydx () const |
Get a list of vectors for the corresponding derivatives. | |
Control numerical methods | |
void | set_mroot (mroot< mm_funct, ubvector, jac_funct > &s_mrp) |
Set solver. | |
void | set_min (min_base<> &s_mp) |
Set minimizer. | |
void | set_stepper (astep_base< ubvector, ubvector, ubvector, ode_funct > &sap) |
Set the adaptive stepper. | |
Public Attributes | |
size_t | buffer_size |
Size of the ODE solution buffer (default ![]() | |
size_t | max_table_size |
Maximum number of lines in the output table (default 400) | |
double | pmax_default |
Default value of maximum pressure for maximum mass star in ![]() ![]() | |
Basic properties | |
double | mass |
Gravitational mass (in ![]() | |
double | rad |
Radius (in km) | |
double | bmass |
Baryonic mass (when computed; in ![]() | |
double | gpot |
Gravitational potential (when computed; unitless) | |
double | last_rjw |
The value of ![]() | |
double | last_f |
The value of ![]() | |
double | domega_rat |
The value of ![]() | |
double | pcent_max |
Maximum value for central pressure in ![]() ![]() | |
Solution parameters | |
These parameters can be changed at any time. | |
bool | reformat_results |
If true, then fixed() and max() reformat results into a o2scl::table_units object. More... | |
double | baryon_mass |
The mass of one baryon. More... | |
bool | ang_vel |
If true, solve for the angular velocity (default false) | |
bool | gen_rel |
Apply general relativistic corrections (default true) | |
bool | calc_gpot |
calculate the gravitational potential (default false) | |
double | step_min |
smallest allowed radial stepsize in km (default 1.0e-4) | |
double | step_max |
largest allowed radial stepsize in km (default 0.05) | |
double | step_start |
initial radial stepsize in km (default 4.0e-3) | |
int | verbose |
control for output (default 1) | |
size_t | max_integ_steps |
Maximum number of integration steps (default 100000) | |
bool | err_nonconv |
If true, call the error handler if the solution does not converge (default true) | |
Mass versus radius parameters | |
double | prbegin |
Beginning pressure in ![]() | |
double | prend |
Ending pressure in ![]() | |
double | princ |
Increment factor for pressure (default 1.1) | |
std::vector< double > | pr_list |
List of pressures at which more information should be recorded. More... | |
Fixed mass parameter | |
double | fixed_pr_guess |
Guess for central pressure in ![]() ![]() | |
Maximum mass profile parameters | |
double | max_begin |
Beginning pressure for maximum mass guess in ![]() | |
double | max_end |
Ending pressure for maximum mass guess in ![]() | |
double | max_inc |
Increment for pressure for maximum mass guess (unitless; default 1.3) | |
Default numerical objects | |
min_brent_gsl | def_min |
The default minimizer. | |
mroot_hybrids< mm_funct, ubvector, ubmatrix, jac_funct > | def_solver |
The default solver. | |
astep_gsl< ubvector, ubvector, ubvector, ode_funct > | def_stepper |
The default adaptive stepper. | |
Protected Member Functions | |
void | column_setup (bool mvsr_mode=false) |
Set up column names and units. More... | |
void | make_unique_name (std::string &col, std::vector< std::string > &cnames) |
Ensure col does not match strings in cnames . More... | |
virtual int | derivs (double x, size_t nv, const ubvector &y, ubvector &dydx) |
The ODE step function. | |
virtual double | max_fun (double maxx) |
The minimizer function to compute the maximum mass. | |
virtual int | integ_star (size_t ndvar, const ubvector &ndx, ubvector &ndy) |
The solver function to compute the stellar profile. | |
Protected Attributes | |
ode_funct | ofm |
ODE function object. | |
interp< ubvector > | iop |
Interpolation object for listed radii in mvsr() | |
bool | integ_star_final |
If true, integ_star() computes all the profile info, otherwise it only computes the gravitational mass. | |
size_t | ix_last |
The last row index in rky. | |
double | schwarz_km |
The schwarzchild radius in km. | |
double | tmass |
Target mass for integ_star() More... | |
double | min_log_pres |
Smallest allowed pressure for integration (default: -100) More... | |
std::shared_ptr< table_units<> > | out_table |
The output table. | |
User EOS | |
eos_tov * | te |
The EOS. | |
bool | eos_set |
True if the EOS has been set. | |
Units for output table | |
std::string | eunits |
Units for energy density. | |
std::string | punits |
Units for pressure. | |
std::string | nunits |
Units for baryon density. | |
double | efactor |
unit conversion factor for energy density (default 1.0) | |
double | pfactor |
unit conversion factor for pressure (default 1.0) | |
double | nfactor |
unit conversion factor for baryon density (default 1.0) | |
Integration storage | |
ubvector | rkx |
Radial coordinate (in kilometers) | |
std::vector< ubvector > | rky |
ODE functions. More... | |
std::vector< ubvector > | rkdydx |
The derivatives of the ODE functions. | |
Numerical methods | |
mroot< mm_funct, ubvector, jac_funct > * | mroot_ptr |
The solver. | |
min_base * | min_ptr |
The minimizer. | |
astep_base< ubvector, ubvector, ubvector, ode_funct > * | as_ptr |
The adaptive stepper. | |
|
protected |
When this function is used by mvsr(), mvsr_mode
is set to true.
|
virtual |
If the target mass is negative, it is interpreted as subtracting from the maximum mass configuration. For a 2.15 solar mass neutron star, target_mass=-0.2
corresponds to 1.95 solar mass configuration.
The variable pmax
is the maximum allowable central pressure in , i.e. the central pressure of the maximum mass star. This ensures that the function does not unintentionally select a configuration on an unstable branch. If
pmax
is greater than or equal to the default value (pmax_default), then the maximum mass star will be explicitly computed first.
|
virtual |
This function constructs a table_units object from the information in rkx, rky, and rkdydx . Note that the table is always constructed by default so this function need not be called unless tov_solve::reformat_results is false
.>
|
protected |
Underscores are added to col
until it matches none of the strings in cnames
.
|
inline |
This function immediately adds four constants to the table, schwarz, Msun, pi
and mproton
.
Definition at line 673 of file tov_solve.h.
double o2scl::tov_solve::baryon_mass |
The mass of one baryon in kg for the total baryon mass calculation (defaults to the proton mass).
Definition at line 522 of file tov_solve.h.
double o2scl::tov_solve::fixed_pr_guess |
This guess is used in the functions fixed().
Definition at line 585 of file tov_solve.h.
|
protected |
This quantity can't be much smaller than -100 since we need to compute numbers near
Definition at line 415 of file tov_solve.h.
double o2scl::tov_solve::pcent_max |
This variable is set by the mvsr()
and max()
functions and used in integ_star() .
Definition at line 502 of file tov_solve.h.
std::vector<double> o2scl::tov_solve::pr_list |
If pressures (in the user-specified units) are added to this vector, then in mvsr(), the radial location, enclosed gravitational mass, and (if o2scl::eos_tov::baryon_column is true) enclosed baryon mass are stored in the table for each central pressure. The associated columns are named r0, gm0, bm0, r1, gm1, bm1,
etc.
Definition at line 574 of file tov_solve.h.
bool o2scl::tov_solve::reformat_results |
Note that mvsr() always places its results in the output table independent of the value of this variable.
Definition at line 516 of file tov_solve.h.
|
protected |
If rky
is viewed as a matrix, then the first column of each row is the gravitational mass in solar masses, and the second column is the natural logarithm of the pressure in . When calc_gpot is true, the next column is the gravitational potential (which is unitless), and when eos_tov::baryon_column is true, the next column is the baryonic mass in
.
Definition at line 432 of file tov_solve.h.
|
protected |
Has a value of zero, unless set to a non-zero value by fixed().
Definition at line 373 of file tov_solve.h.
Documentation generated with Doxygen. Provided under the
GNU Free Documentation License (see License Information).