Public Types | Public Attributes | Protected Member Functions | Protected Attributes | List of all members
o2scl::tov_solve Class Reference

Class to solve the Tolman-Oppenheimer-Volkov equations. More...

#include <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< ubmatrixubmatrix_column
typedef boost::numeric::ublas::matrix_row< ubmatrixubmatrix_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, double pmax=1.0e20)
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 ubvectorget_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_funct11, ubvector, jac_funct11 > &s_mrp)
 Set solver.
void set_min (min_base<> &s_mp)
 Set minimizer.
void set_stepper (astep_base< ubvector, ubvector, ubvector, ode_funct11 > &sap)
 Set the adaptive stepper.

Public Attributes

size_t buffer_size
 Size of the ODE solution buffer (default $ 10^5 $)
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 $ \mathrm{M}_{\odot}/\mathrm{km}^3 $.
Basic properties
double mass
 Gravitational mass.
double rad
double bmass
 Baryonic mass (when computed)
double gpot
 Gravitational potential (when computed)
double last_rjw
 The value of $ r^4 j df / dr $.
double last_f
 The value of $ \bar{\omega} / \Omega $.
double domega_rat
 The value of $ d (\bar{\omega}/\Omega)/dr $ at the surface (when computed)
double pcent_max
 Maximum value for central pressure in $ \mathrm{M}_{\odot}/\mathrm{km}^3 $ (default $ 10^{20} $ ) More...
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 and the enclosed baryon mass (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 $ \mathrm{M}_{\odot}/\mathrm{km}^3 $ (default 7.0e-7)
double prend
 Ending pressure in $ \mathrm{M}_{\odot}/\mathrm{km}^3 $ (default 8.0e-3)
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 $ \mathrm{M}_{\odot}/\mathrm{km}^3 $ (default $ 5.2 \times 10^{-5} $) More...
Maximum mass profile parameters
double max_begin
 Beginning pressure for maximum mass guess (default 7.0e-5)
double max_end
 Ending pressure for maximum mass guess (default 5.0e-3)
double max_inc
 Increment for pressure for maximum mass guess (default 1.3)
Default numerical classes
min_brent_gsl def_min
 The default minimizer.
mroot_hybrids< mm_funct11, ubvector, ubmatrix, jac_funct11def_solver
 The default solver.
astep_gsl< ubvector, ubvector, ubvector, ode_funct11def_stepper
 The default adaptive stepper.

Static Public Attributes

Convergence information flags
static const int fixed_solver_failed =128
static const int fixed_integ_star_failed =256
static const int fixed_wrong_mass =512
static const int max_minimizer_failed =1024
static const int max_integ_star_failed =2048
static const int mvsr_integ_star_failed =4096
static const int ang_vel_failed =8192
static const int cent_press_large =16384
static const int cent_press_neg =32768
static const int over_max_steps =65536
static const int last_step_large =131072

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_funct11 ofm
 ODE function object.
interp< ubvectoriop
 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
 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< ubvectorrky
 ODE functions. More...
std::vector< ubvectorrkdydx
 The derivatives of the ODE functions.
Numerical methods
mroot< mm_funct11, ubvector, jac_funct11 > * mroot_ptr
 The solver.
 The minimizer.
astep_base< ubvector, ubvector, ubvector, ode_funct11 > * as_ptr
 The adaptive stepper.

Detailed Description

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(), or tov_solve::max() to compute the mass versus radius curve, the profile of a star of a given fixed mass, or the profile of the maximum mass star. These functions all generate results in the form of a table_units object, which can be obtained from tov_solve::get_results().

Screen output:

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. The intial guess to the solver is always the value of fixed_pr_guess, which defaults to $ 5.2 \times 10^{-5}~\mathrm{Msun}/\mathrm{km}^3 $ . Alternatively, the user can specify the central pressure of the maximum mass star so that it does not have to be computed.

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).

Output tables

The functions tov_solve::fixed() and tov_solve::max() produce output tables which represent the profile of the neutron star of the requested mass. The columns are

By default, this class operates with energy density and pressure in units of $ \mathrm{M}_{\odot}/\mathrm{km}^3 $ and baryon density in units of $ \mathrm{fm}^{-3} $.

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:

The list of hard-coded units for baryon density are:

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.


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 $ 10^8 $.


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 $ \bar{\omega} $ (see the section in the User's Guide called Solution of the Tolman-Oppenheimer-Volkov equations ) is independent of the relative scale for $ \bar{\omega} $ and $ j $ . First, one rescales $ \bar{\omega} $ and rewrites everything in terms of $ f\equiv \bar{\omega}/\Omega $ and $ g \equiv r^4 j~df/dr $ . Then, pick a central pressure, $ m(r=0)=g(r=0)=0 $, arbitrary values for $ \Phi(r=0) $ and $ f(r=0) $, and integrate

\begin{eqnarray*} \frac{d P}{d r} &=& - \frac{G \varepsilon m}{r^2} \left( 1+\frac{P}{\varepsilon}\right) \left( 1+\frac{4 \pi P r^3}{m} \right) \left( 1-\frac{2 G m}{r}\right)^{-1} \nonumber \\ \frac{d m}{d r} &=& 4 \pi r^2 \varepsilon \nonumber \\ \frac{d \Phi}{d r} &=& - \frac{1}{\varepsilon} \frac{ d P}{d r} \left(1+\frac{P}{\varepsilon}\right)^{-1} \nonumber \\ \frac{d g}{dr} &=& -4 r^3 \frac{d j}{dr} f \nonumber \\ \frac{d f}{dr} &=& \frac{g}{r^4 j} \end{eqnarray*}

Afterwards, shift $ \Phi $ by a constant to ensure the correct value at $ r=R $, and multiply $ g $ by a constant to ensure that $ g=r^4 j (df/dr) $ holds for the new potential $ \Phi $. Then, multiply $ f $ by a constant to ensure that

\[ f(r=R) + \frac{R}{3} \left(\frac{df}{dr}\right)_{r=R} = 1 \]

The functions $ f $ and $ g $ 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),

\[ I = \frac{R^4}{6 G} \left.\frac{df}{dr}\right|_{r=R} \]

where the last fraction is stored in domega_rat .

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.

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.

The function o2scl::tov_solve::integ_star() returns 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 265 of file tov_solve.h.

Member Function Documentation

◆ column_setup()

void o2scl::tov_solve::column_setup ( bool  mvsr_mode = false)

When this function is used by mvsr(), mvsr_mode is set to true.

◆ fixed()

virtual int o2scl::tov_solve::fixed ( double  target_mass,
double  pmax = 1.0e20 

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 $ \mathrm{M}_{\odot}/\mathrm{km}^3 $, 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.

◆ make_table()

virtual void o2scl::tov_solve::make_table ( )

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.

◆ make_unique_name()

void o2scl::tov_solve::make_unique_name ( std::string &  col,
std::vector< std::string > &  cnames 

Underscores are added to col until it matches none of the strings in cnames.

◆ set_table()

void o2scl::tov_solve::set_table ( std::shared_ptr< table_units<> >  t)

This function immediately adds four constants to the table, schwarz, Msun, pi and mproton.

Definition at line 597 of file tov_solve.h.

Member Data Documentation

◆ baryon_mass

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 455 of file tov_solve.h.

◆ fixed_pr_guess

double o2scl::tov_solve::fixed_pr_guess

This guess is used in the functions fixed().

Definition at line 518 of file tov_solve.h.

◆ min_log_pres

double o2scl::tov_solve::min_log_pres

This quantity can't be much smaller than -100 since we need to compute numbers near $ e^{-\mathrm{min\_log\_pres}} $

Idea for Future:
Replace this with the proper value from std::limits?

Definition at line 348 of file tov_solve.h.

◆ pcent_max

double o2scl::tov_solve::pcent_max

This variable is set by the mvsr() and max() functions and used in integ_star() .

Definition at line 435 of file tov_solve.h.

◆ pr_list

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 507 of file tov_solve.h.

◆ reformat_results

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 449 of file tov_solve.h.

◆ rky

std::vector<ubvector> o2scl::tov_solve::rky

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 $ \mathrm{M}_{\odot}/km^3 $ . 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 $ \mathrm{M}_{\odot} $.

Definition at line 365 of file tov_solve.h.

◆ tmass

double o2scl::tov_solve::tmass

Has a value of zero, unless set to a non-zero value by fixed().

Definition at line 306 of file tov_solve.h.

The documentation for this class was generated from the following file:

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).