eos_nse.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file eos_nse.h
24  \brief File defining \ref o2scl::eos_nse
25 */
26 #ifndef O2SCL_NSE_EOS_H
27 #define O2SCL_NSE_EOS_H
28 
29 #include <o2scl/classical.h>
30 #include <o2scl/constants.h>
31 #include <o2scl/nucdist.h>
32 #include <o2scl/mm_funct.h>
33 #include <o2scl/mroot_hybrids.h>
34 #include <o2scl/mmin_simp2.h>
35 
36 #ifndef DOXYGEN_NO_O2NS
37 namespace o2scl {
38 #endif
39 
40  /** \brief Equation of state for nuclei in statistical equilibrium
41 
42  This class computes the composition of matter in nuclear
43  statistical equilibrium. The chemical potential of a nucleus X
44  with proton number \f$ Z_X \f$ and neutron number \f$ N_X \f$ is
45  given by
46  \f[
47  \mu_X = N \mu_n + Z \mu_p - E_{\mathrm{bind},X}
48  \f]
49  where \f$ \mu_n \f$ and \f$ \mu_p \f$ are the neutron and proton
50  chemical potentials and \f$ E_{\mathrm{bind},X} \f$ is the
51  binding energy of the nucleus. The chemical potentials are
52  assumed to be in units of \f$ \mathrm{fm}^{-1} \f$.
53 
54  The baryon number density and electron fraction are then given
55  by
56  \f[
57  n_B = \sum_X n_{X} (N_X + Z_X) \qquad Y_e n_B = \sum_X n_X Z_X
58  \f]
59  where \f$ n_X \f$ is the number density which is determined from
60  the chemical potential above.
61 
62  The nuclei in specified in the parameter named \c nd, must have
63  their proton number, neutron number, mass number, binding
64  energy, and spin degeracy already specified. This class
65  implicitly assumes that the nuclei are non-interacting and that
66  the values of <tt>part::inc_rest_mass</tt> are false. The
67  chemical potential arguments also do not include the rest mass.
68  The nuclear rest mass is presumed to be \f$ Z_X m_p + N_X m_n
69  \f$.
70 
71  The function \ref calc_density() attempts to solve for the
72  neutron and proton chemical potentials given the neutron and
73  proton densities. However, this is relatively difficult. At low
74  enough temperatures, \f$ n(\mu) \f$ is a staircase-like function
75  with alernating regions which are very flat and or nearly
76  vertical. For this reason, derivative-based methods often fail
77  without extremely good guesses. The current method of solution
78  combines \ref make_guess(), \ref density_min() and \ref
79  direct_solve() in order to obtain the solution.
80 
81  Note also that \ref calc_density() will fail if there are
82  no nuclei in the distribution which equal, or surround the
83  requested value of \f$ Y_e=n_p/(n_n+n_p) \f$ determined from \c nn and
84  \c np .
85  */
86  class eos_nse {
87 
88  public:
89 
91 
92 #ifndef DOXYGEN_INTERNAL
93 
94  protected:
95 
96  /// Function to solve to match neutron and proton densities
97  int solve_fun(size_t nv, const ubvector &x, ubvector &y,
98  double nn, double np, double T,
99  std::vector<o2scl::nucleus> &nd);
100 
101  /// Function to minimize to match neutron and proton densities
102  double minimize_fun(size_t nv, const ubvector &x, double T,
103  double nn, double np, o2scl::thermo &th,
104  std::vector<o2scl::nucleus> &nd);
105 
106  /// Solver
108 
109  /// Minimizer
111 
112  /// Compute particle properties assuming classical thermodynamics
114 
115 #endif
116 
117  public:
118 
119  eos_nse();
120 
121  /// \name Basic usage
122  //@{
123  /// Verbosity parameter (default 1)
124  int verbose;
125 
126  /** \brief If true, call the error handler if calc_density() does
127  not converge (default true)
128  */
130 
131  /** \brief Calculate the equation of state as a function of the
132  chemical potentials
133 
134  Given \c mun, \c mup and \c T, this computes the composition
135  (the individual densities are stored in the distribution \c
136  nd), the neutron number density \c nn, and the proton number
137  density \c np. Note that the densities can be infinite if
138  the chemical potentials are sufficiently large.
139 
140  This function does not use the solver or the minimizer.
141  */
142  void calc_mu(double mun, double mup, double T, double &nn,
143  double &np, thermo &th, std::vector<nucleus> &nd);
144 
145  /** \brief Calculate the equation of state as a function of the densities
146 
147  Given the neutron number density \c nn in \f$ \mathrm{fm}^{-3}
148  \f$, the proton number density \c np and the temperature \c T
149  in \f$ \mathrm{fm}^{-1} \f$, this computes the composition
150  (the individual densities are stored in the distribution \c
151  nd) and the chemical potentials are given in \c mun and \c mup
152  . The nuclei in \c nd must have their proton number, neutron
153  number, atomic number, binding energy, and spin degeracy
154  already specified.
155 
156  This function uses \ref make_guess(), \ref direct_solve(),
157  and \ref density_min(), to self-consistently compute the
158  chemical potentials.
159  */
160  int calc_density(double nn, double np, double T, double &mun,
161  double &mup, thermo &th, std::vector<nucleus> &nd);
162  //@}
163 
164  /// \name Tools for fixing chemical potentials from the densities
165  //@{
166  /// The maximum number of iterations for \ref make_guess() (default 60)
168 
169  /** \brief The initial stepsize for the chemical potentials relative
170  to the temperature (default \f$ 10^5 \f$ )
171  */
173 
174  /** \brief Find values for the chemical potentials which ensure
175  that the densities are within a fixed range
176 
177  This function improves initial guesses for the chemical
178  potentials in order to ensure the densities are within a
179  specified range. It can sometimes even succeed when the
180  chemical potentials are so far off as to make the densities
181  infinite or zero. This function is used by \ref calc_density()
182  to improve the initial guesses for the chemical potentials if
183  necessary.
184 
185  The algorithm can fail in several different ways. This is
186  particularly likely if the density range specified by \c
187  nn_min, \c nn_max, \c np_min, and \c np_max is small. This
188  function ignores the value of \ref err_nonconv, and throws an
189  exception on failure only if \c err_on_fail is true (which is
190  the default).
191  */
192  int make_guess(double &mun, double &mup, double T,
193  thermo &th, std::vector<nucleus> &nd,
194  double nn_min=1.0e-20, double nn_max=1.0e8,
195  double np_min=1.0e-20, double np_max=1.0e8,
196  bool err_on_fail=true);
197 
198  /** \brief Obtain chemical potentials from densities directly
199  using a solver
200 
201  This function often requires extremely good guesses for the
202  chemical potentials, especially at low temperatures.
203  */
204  int direct_solve(double nn, double np, double T,
205  double &mun, double &mup, thermo &th,
206  std::vector<nucleus> &nd);
207 
208  /** \brief Obtain chemical potentials from densities
209  using a minimizer
210 
211  This function often requires extremely good guesses for the
212  chemical potentials, especially at low temperatures. By
213  default, this calls the minimizer five times, as this seems to
214  improve convergence using the default minimizer. By default,
215  the value of \ref o2scl::mmin_base::err_nonconv is set to
216  false for \ref def_mmin .
217  */
218  int density_min(double nn, double np, double T,
219  double &mun, double &mup, thermo &th,
220  std::vector<nucleus> &nd);
221  //@}
222 
223  /// \name Numerical methods
224  //@{
225  /// Default solver
227 
228  /// Default minimizer
230 
231  /** \brief Set the solver for use in \ref direct_solve()
232  */
233  void set_mroot(mroot<> &rp) {
234  mroot_ptr=&rp;
235  return;
236  }
237 
238  /** \brief Set the minimizer for use in \ref density_min()
239  */
240  void set_mmin(mmin_base<> &mp) {
241  mmin_ptr=&mp;
242  return;
243  }
244  //@}
245 
246  };
247 
248 #ifndef DOXYGEN_NO_O2NS
249 }
250 #endif
251 
252 #endif
void set_mmin(mmin_base<> &mp)
Set the minimizer for use in density_min()
Definition: eos_nse.h:240
size_t make_guess_iters
The maximum number of iterations for make_guess() (default 60)
Definition: eos_nse.h:167
bool err_nonconv
If true, call the error handler if calc_density() does not converge (default true) ...
Definition: eos_nse.h:129
mmin_base * mmin_ptr
Minimizer.
Definition: eos_nse.h:110
int make_guess(double &mun, double &mup, double T, thermo &th, std::vector< nucleus > &nd, double nn_min=1.0e-20, double nn_max=1.0e8, double np_min=1.0e-20, double np_max=1.0e8, bool err_on_fail=true)
Find values for the chemical potentials which ensure that the densities are within a fixed range...
double make_guess_init_step
The initial stepsize for the chemical potentials relative to the temperature (default ) ...
Definition: eos_nse.h:172
int direct_solve(double nn, double np, double T, double &mun, double &mup, thermo &th, std::vector< nucleus > &nd)
Obtain chemical potentials from densities directly using a solver.
mmin_simp2 def_mmin
Default minimizer.
Definition: eos_nse.h:229
void set_mroot(mroot<> &rp)
Set the solver for use in direct_solve()
Definition: eos_nse.h:233
int density_min(double nn, double np, double T, double &mun, double &mup, thermo &th, std::vector< nucleus > &nd)
Obtain chemical potentials from densities using a minimizer.
classical cla
Compute particle properties assuming classical thermodynamics.
Definition: eos_nse.h:113
int calc_density(double nn, double np, double T, double &mun, double &mup, thermo &th, std::vector< nucleus > &nd)
Calculate the equation of state as a function of the densities.
double minimize_fun(size_t nv, const ubvector &x, double T, double nn, double np, o2scl::thermo &th, std::vector< o2scl::nucleus > &nd)
Function to minimize to match neutron and proton densities.
int verbose
Verbosity parameter (default 1)
Definition: eos_nse.h:124
mroot * mroot_ptr
Solver.
Definition: eos_nse.h:107
mroot_hybrids def_mroot
Default solver.
Definition: eos_nse.h:226
Equation of state for nuclei in statistical equilibrium.
Definition: eos_nse.h:86
void calc_mu(double mun, double mup, double T, double &nn, double &np, thermo &th, std::vector< nucleus > &nd)
Calculate the equation of state as a function of the chemical potentials.
int solve_fun(size_t nv, const ubvector &x, ubvector &y, double nn, double np, double T, std::vector< o2scl::nucleus > &nd)
Function to solve to match neutron and proton densities.

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