eos_nse_full.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2014-2020, 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_full.h
24  \brief File defining \ref o2scl::eos_nse_full
25 */
26 #ifndef EOS_NSE_FULL_H
27 #define EOS_NSE_FULL_H
28 
29 #include <boost/numeric/ublas/vector.hpp>
30 
31 #include <o2scl/constants.h>
32 
33 #include <o2scl/classical.h>
34 #include <o2scl/fermion_rel.h>
35 #include <o2scl/fermion_deriv_rel.h>
36 
37 #include <o2scl/nucmass_densmat.h>
38 #include <o2scl/mroot_hybrids.h>
39 #include <o2scl/root_cern.h>
40 
41 #include <o2scl/eos_had_skyrme.h>
42 
43 #include <o2scl/mmin_simp2.h>
44 
45 #ifndef DOXYGEN_NO_O2NS
46 namespace o2scl {
47 #endif
48 
49  /** \brief EOS for nuclear statistical equilibrium with interactions
50 
51  This class is experimental.
52 
53  For the verbose parameter, generally 0 means no output, 1 means
54  the function will output the composition and thermodynamics for
55  the first 10 or so nuclei in the distribution, and 2 means the
56  function will output the entire distribution.
57 
58  This class retains the usual mechanism using \ref err_nonconv to
59  handle what to do if one of the functions does not converge. In
60  addition, \ref calc_density_fixnp() and \ref
61  calc_density_noneq() return \ref invalid_config for invalid
62  configurations, which sometimes occur during normal execution.
63  Since these invalid configurations are 'normal', they do not
64  cause the error handler to be called, independent of the value
65  of \ref err_nonconv . Practically, this means the end-user must
66  check the return value of these two functions every time they
67  are called.
68 
69  This class presumes that electrons include their rest mass,
70  but nucleons and nuclei do not. The error handler is called
71  by some functions if this is not the case (determined
72  by the values in <tt>o2scl::part::inc_rest_mass</tt>).
73 
74  \todo I don't think inc_lept_phot=false works because then
75  all WS cells have infinite size because of no electrons.
76  For the moment, this variable is protected to discourage
77  the user from changing it.
78 
79  \future There is a bit of duplication between calc_density_noneq()
80  and calc_density_fixnp() which could be streamlined.
81  \future Add fermion and boson statistics to the nuclei in the
82  distribution.
83  */
84  class eos_nse_full {
85 
86  public:
87 
90 
91  protected:
92 
93  /** \brief Check the \ref o2scl::dense_matter object to
94  see if the rest masses are correctly included or not,
95  etc.
96  */
97  virtual void check_dm(o2scl::dense_matter &dm);
98 
99  /** \brief Output a \ref o2scl::dense_matter object
100  according to the setting of \ref verbose
101  for function specified in \c func_name .
102  */
103  virtual void verb_output(o2scl::dense_matter &dm,
104  std::string func_name);
105 
106  /** \brief If true, include electrons and photons (default true)
107  */
109 
110  /// Compute particle properties assuming classical thermodynamics
111  o2scl::classical_thermo cla;
112 
113  /// Relativistic fermions with derivatives
114  o2scl::fermion_deriv_rel snf;
115 
116  /// Mass formula (points to \ref nuc_dens by default)
118 
119  /** \brief The full distribution of all nuclei to consider
120 
121  \note Currently, the \c ad variable doesn't do much, but
122  it's important to leave this in as future functions may
123  want to automatically adjust the distribution
124  */
125  std::vector<o2scl::nucleus> *ad;
126 
127  /** \brief Solve for charge neutrality assuming the specified
128  electron chemical potential and proton number density.
129 
130  This function returns
131  \f[
132  \left[n_p-n_e(\mu_e)-n_{\mu}(\mu_{\mu}=\mu_e)\right]/n_p
133  \f]
134  using \ref relf.
135  */
136  virtual double charge_neutrality(double mu_e, double np_tot,
137  dense_matter &dm);
138 
139  /** \brief Compute the free energy from a vector of densities
140  of the nuclei
141 
142  This calls \ref calc_density_noneq() and then returns the free
143  energy. The vector \c n_nuc and the distribution \c dm.dist
144  must both have the same size. The nuclear densities are
145  taken from \c n_nuc and the proton and neutron densities
146  are determined automatically from subtracting the density
147  contributions of nuclei from the total neutron and proton
148  densities as determined in \ref o2scl::dense_matter::nB
149  and \ref o2scl::dense_matter::Ye .
150 
151  If the call to \ref calc_density_noneq() returns a non-zero
152  value, e.g. because of an invalid configuration,
153  then the value \f$ 10^{4} \f$ is returned.
154  */
155  virtual double free_energy(const ubvector &n_nuc, dense_matter &dm);
156 
157  /// Nucleonic EOS (0 by default)
159 
160  public:
161 
162  eos_nse_full();
163 
164  /** \brief The integer return value which indicates an invalid
165  configuration
166  */
167  static const int invalid_config=-10;
168 
169  /// \name Various settings
170  //@{
171  /// Verbose parameter
172  int verbose;
173 
174  /** \brief If true, call the error handler if calc_density() does
175  not converge (default true)
176  */
178 
179  /** \brief If true, include dripped protons and
180  neutrons in the nuclear mass (default true)
181  */
183 
184  /// If true, include muons (default false)
186  //@}
187 
188  /** \brief Function which is solved by \ref calc_density_saha()
189 
190  This function takes two inputs, the neutron and proton
191  densities, and solves to ensure that \ref
192  dense_matter::baryon_density() matches \ref
193  o2scl::dense_matter::nB and \ref
194  o2scl::dense_matter::electron_fraction() matches \ref
195  dense_matter::Ye.
196 
197  This function calls \ref calc_density_fixnp() .
198  */
199  virtual int solve_fixnp(size_t n, const ubvector &x, ubvector &y,
200  dense_matter &dm, bool from_densities=true);
201 
202  /** \brief Solve matter at fixed chemical potential by
203  bracketing
204  */
205  virtual int bracket_mu_solve(double &mun_low, double &mun_high,
206  double &mup_low, double &mup_high,
207  dense_matter &dm);
208 
209  /** \brief Fix electron fraction by varying proton chemical
210  potential
211 
212  At some fixed values of <tt>dm.Ye</tt> and <tt>dm.nB</tt>,
213  given a value of \f$ \mu_p \f$, and given an initial bracket
214  for \f$ \mu_n \f$ (stored in <tt>mun_low</tt> and
215  <tt>mun_high</tt>), this function attempts to find the value
216  of \f$ \mu_n \f$ which ensures that the baryon density in
217  nuclei matches that in <tt>dm.nB</tt> using a bracketing
218  solver. It then returns the difference between the value of
219  the proton fraction in nuclei and the value in <tt>dm.Ye</tt>.
220 
221  If <tt>mun_low</tt> and <tt>mun_high</tt> do not bracket the
222  correct value of \f$ \mu_n \f$, this function attempts to
223  modify them to give a proper bracket for the root. The
224  finaly value of \f$ \mu_n \f$ is stored in <tt>dm.n.mu</tt>.
225 
226  Currently, the values of <tt>dm.n.n</tt> and <tt>dm.p.n</tt>
227  are ignored and set to zero.
228  */
229  double mup_for_Ye(double mup, double &mun_low,
230  double &mun_high, dense_matter &dm);
231 
232  /** \brief Fix the baryon density by varying the neutron
233  chemical potential
234 
235  Given a value of \f$ \mu_n \f$ (the value in <tt>dm.n.mu</tt>
236  is ignored), this function computes the baryon density
237  in nuclei and returns the difference between this value
238  and that stored in <tt>dm.nB</tt>.
239 
240  Currently, the values of <tt>dm.n.n</tt> and <tt>dm.p.n</tt>
241  are ignored and set to zero.
242  */
243  virtual double solve_mun(double mun, dense_matter &dm);
244 
245  /** \brief Compute the properties of matter from the densities,
246  not presuming equilibrium
247 
248  The values of <tt>dm.nB</tt> and <tt>dm.Ye</tt> are ignored
249  and unchanged by this function. The electron and muon density
250  are determined by charged neutrality and assuming their
251  chemical potentials are equal. Photons are always included.
252 
253  If the nuclear densities are all zero, then this just
254  returns nuclear matter with leptons and photons as
255  determined by charge neutrality.
256 
257  This function is designed to return non-zero values for
258  invalid configurations and can return the value
259  \ref invalid_config without calling the error handler,
260  independent of the value of \ref err_nonconv .
261 
262  Possible invalid configurations are:
263  - negative nucleon or nucleus densities, or
264  - proton radii larger than WS cell radii, i.e.
265  \f$ (0.08 - n_p) / (n_e+n_{\mu}-n_p) < 1 \f$ or
266  \f$ n_p > 0.08 \f$ .
267  */
268  virtual int calc_density_noneq(dense_matter &dm);
269 
270  /** \brief Compute the properties of matter from
271  neutron and proton densities, using the Saha equation
272 
273  If the parameter <tt>from_densities</tt> is true, then this
274  computes nucleonic matter using the neutron and proton
275  densities stored in <tt>dm.n.n</tt> and <tt>dm.p.n</tt>.
276  Otherwise, nucleonic matter is computed using the chemical
277  potential stored in <tt>dm.n.mu</tt> and <tt>dm.p.mu</tt>.
278  Either way, electrons are computed assuming their density is
279  given from \ref o2scl::dense_matter::nB and \ref
280  o2scl::dense_matter::Ye. Muons are added assuming their
281  chemical potential is equal to the electron chemical
282  potential. Finally, the Saha equation is used to determine the
283  nuclear chemical potentials and this gives the nuclear
284  densities.
285 
286  This function only works when \ref inc_prot_coul is
287  <tt>false</tt>.
288 
289  The values in \ref o2scl::dense_matter::nB and \ref
290  o2scl::dense_matter::Ye are unchanged by this function. Note
291  that, after this function completes, the value returned by
292  \ref o2scl::dense_matter::baryon_density() will not
293  necessarily be the same as that stored in \ref
294  o2scl::dense_matter::nB (and similarly for the electron
295  fraction).
296 
297  This function is designed to return non-zero values for
298  invalid configurations and can return the value
299  \ref invalid_config without calling the error handler,
300  independent of the value of \ref err_nonconv .
301 
302  Possible invalid configurations are:
303  - negative nucleon densities, or
304  - proton radii larger than WS cell radii, i.e.
305  \f$ (0.08 - n_p) / (n_e+n_{\mu}-n_p) < 1 \f$ or
306  \f$ n_p > 0.08 \f$ .
307  */
308  virtual int calc_density_fixnp(dense_matter &dm, bool from_densities=true);
309 
310  /** \brief Compute the free energy for a fixed composition
311  by minimization
312 
313  Given a fixed baryon density (dm.nB), electron fraction
314  (dm.Ye), temperature (dm.T), this minimizes the free energy
315  over the densities of the nuclei currently present in the
316  distribution. The neutron and proton drip densities are
317  determined by ensuring that the baryon density and electron
318  fraction are correctly reproduced. The function which
319  is minimized is \ref free_energy() .
320 
321  \note This function currently only performs a very simple
322  minimization and currently works in only limited
323  circumstances.
324  */
325  virtual int calc_density_by_min(dense_matter &dm);
326 
327  /** \brief Compute properties of matter for baryon density and
328  electron fraction using the Saha equation
329 
330  This function solves the function specified by \ref
331  solve_fixnp() using the current values of <tt>dm.n.n</tt> and
332  <tt>dm.p.n</tt> as initial guesses.
333  */
334  virtual int calc_density_saha(dense_matter &dm);
335 
336  /** \brief Output properties of a \ref o2scl::dense_matter object to
337  std::cout
338 
339  This function was particularly designed for comparing results
340  with \ref o2scl::eos_sn_base derived classes.
341 
342  If output level is 0, then just the basic quantities are
343  output without any information about the distribution. If
344  output_level is 1, then only about 10 nuclei in the
345  distribution are output, and if output_level is 2,
346  then all nuclei in the distribution are output.
347  */
348  virtual void output(dense_matter &dm, int output_level);
349 
350  /** \brief Adjust the particle densities to match specified
351  density and electron fraction
352 
353  This function attempts to match the nuclear and nucleon
354  densities so that the baryon density and electron fraction are
355  equal to those specified in \ref o2scl::dense_matter::nB and
356  \ref o2scl::dense_matter::Ye .
357  */
358  virtual int density_match(dense_matter &dm);
359 
360  /** \brief Relativistic fermions
361 
362  \comment
363  Must currently be public for tcan/ecn.cpp.
364  \endcomment
365  */
366  o2scl::fermion_rel relf;
367 
368  /// \name Nuclei and nuclear masses
369  //@{
370  /// Compute nuclei in dense matter
372 
373  /** \brief Set nuclear mass formula
374  */
376  massp=&m;
377  return;
378  }
379 
380  /** \brief Set distribution of nuclei
381  */
382  void set_dist(std::vector<o2scl::nucleus> &dist) {
383  ad=&dist;
384  return;
385  }
386  //@}
387 
388  /// \name Nucleonic matter EOS
389  //@{
390  /** \brief Set homogeneous matter EOS
391  */
393  ehtp=&e;
394  return;
395  }
396 
397  /** \brief Get homogeneous matter EOS
398 
399  This function calls the error handler if no EOS has been set
400  */
402  if (ehtp==0) {
403  O2SCL_ERR2("Homogeneous matter EOS not specified in ",
404  "eos_nse_full::get_eos().",exc_efailed);
405  }
406  return *ehtp;
407  }
408 
409  /** \brief Return true if an EOS was specified
410  */
411  bool is_eos_set() {
412  if (ehtp==0) return false;
413  return true;
414  }
415  //@}
416 
417  /// \name Numerical methods
418  //@{
419  /// The default minimizer
421 
422  /// Default solver
424 
425  /// Lepton solver
427  //@}
428 
429 #ifdef O2SCL_NEVER_DEFINED
430 
431  /** \brief Desc
432 
433  This was intended to be a version of calc_density_by_min()
434  which optimized the composition, but it doesn't really work
435  yet.
436  */
437  int calc_density(dense_matter &dm);
438 
439 #endif
440 
441  };
442 
443 #ifndef DOXYGEN_NO_O2NS
444 }
445 #endif
446 
447 #endif
o2scl::mmin_simp2
o2scl::eos_nse_full::set_dist
void set_dist(std::vector< o2scl::nucleus > &dist)
Set distribution of nuclei.
Definition: eos_nse_full.h:382
boost::numeric::ublas::matrix< double >
o2scl::eos_nse_full::calc_density_fixnp
virtual int calc_density_fixnp(dense_matter &dm, bool from_densities=true)
Compute the properties of matter from neutron and proton densities, using the Saha equation.
o2scl::eos_nse_full::def_mmin
o2scl::mmin_simp2 def_mmin
The default minimizer.
Definition: eos_nse_full.h:420
o2scl::eos_nse_full::free_energy
virtual double free_energy(const ubvector &n_nuc, dense_matter &dm)
Compute the free energy from a vector of densities of the nuclei.
boost::numeric::ublas::vector< double >
o2scl::eos_nse_full::cla
o2scl::classical_thermo cla
Compute particle properties assuming classical thermodynamics.
Definition: eos_nse_full.h:111
o2scl::eos_nse_full::calc_density_noneq
virtual int calc_density_noneq(dense_matter &dm)
Compute the properties of matter from the densities, not presuming equilibrium.
exc_efailed
exc_efailed
o2scl::eos_nse_full::inc_prot_coul
bool inc_prot_coul
If true, include dripped protons and neutrons in the nuclear mass (default true)
Definition: eos_nse_full.h:182
o2scl::eos_nse_full::charge_neutrality
virtual double charge_neutrality(double mu_e, double np_tot, dense_matter &dm)
Solve for charge neutrality assuming the specified electron chemical potential and proton number dens...
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
o2scl::eos_nse_full::massp
o2scl::nucmass_densmat * massp
Mass formula (points to nuc_dens by default)
Definition: eos_nse_full.h:117
o2scl::eos_nse_full::mup_for_Ye
double mup_for_Ye(double mup, double &mun_low, double &mun_high, dense_matter &dm)
Fix electron fraction by varying proton chemical potential.
o2scl::eos_nse_full::solve_fixnp
virtual int solve_fixnp(size_t n, const ubvector &x, ubvector &y, dense_matter &dm, bool from_densities=true)
Function which is solved by calc_density_saha()
o2scl::eos_nse_full::calc_density_by_min
virtual int calc_density_by_min(dense_matter &dm)
Compute the free energy for a fixed composition by minimization.
o2scl::eos_nse_full::snf
o2scl::fermion_deriv_rel snf
Relativistic fermions with derivatives.
Definition: eos_nse_full.h:114
o2scl::eos_nse_full::ad
std::vector< o2scl::nucleus > * ad
The full distribution of all nuclei to consider.
Definition: eos_nse_full.h:125
o2scl::eos_nse_full::density_match
virtual int density_match(dense_matter &dm)
Adjust the particle densities to match specified density and electron fraction.
o2scl::eos_nse_full::relf
o2scl::fermion_rel relf
Relativistic fermions.
Definition: eos_nse_full.h:366
o2scl::eos_nse_full::check_dm
virtual void check_dm(o2scl::dense_matter &dm)
Check the o2scl::dense_matter object to see if the rest masses are correctly included or not,...
o2scl::eos_nse_full::def_mroot
mroot_hybrids def_mroot
Default solver.
Definition: eos_nse_full.h:423
o2scl::eos_had_temp_base
A finite temperature hadronic EOS [abstract base].
Definition: eos_had_base.h:967
o2scl::eos_nse_full::bracket_mu_solve
virtual int bracket_mu_solve(double &mun_low, double &mun_high, double &mup_low, double &mup_high, dense_matter &dm)
Solve matter at fixed chemical potential by bracketing.
o2scl::eos_nse_full::output
virtual void output(dense_matter &dm, int output_level)
Output properties of a o2scl::dense_matter object to std::cout.
o2scl::eos_nse_full
EOS for nuclear statistical equilibrium with interactions.
Definition: eos_nse_full.h:84
o2scl::eos_nse_full::verb_output
virtual void verb_output(o2scl::dense_matter &dm, std::string func_name)
Output a o2scl::dense_matter object according to the setting of verbose for function specified in fun...
o2scl::eos_nse_full::set_mass
void set_mass(o2scl::nucmass_densmat &m)
Set nuclear mass formula.
Definition: eos_nse_full.h:375
o2scl::mroot_hybrids
o2scl::nucmass_densmat
o2scl::eos_nse_full::def_root
root_cern def_root
Lepton solver.
Definition: eos_nse_full.h:426
o2scl::eos_nse_full::is_eos_set
bool is_eos_set()
Return true if an EOS was specified.
Definition: eos_nse_full.h:411
o2scl::root_cern
o2scl::eos_nse_full::calc_density_saha
virtual int calc_density_saha(dense_matter &dm)
Compute properties of matter for baryon density and electron fraction using the Saha equation.
o2scl::eos_nse_full::invalid_config
static const int invalid_config
The integer return value which indicates an invalid configuration.
Definition: eos_nse_full.h:167
o2scl::eos_nse_full::inc_lept_phot
bool inc_lept_phot
If true, include electrons and photons (default true)
Definition: eos_nse_full.h:108
o2scl::dense_matter
o2scl::eos_nse_full::err_nonconv
bool err_nonconv
If true, call the error handler if calc_density() does not converge (default true)
Definition: eos_nse_full.h:177
o2scl::eos_nse_full::set_eos
void set_eos(o2scl::eos_had_temp_base &e)
Set homogeneous matter EOS.
Definition: eos_nse_full.h:392
o2scl::eos_nse_full::nuc_dens
o2scl::nucmass_densmat nuc_dens
Compute nuclei in dense matter.
Definition: eos_nse_full.h:371
o2scl::eos_nse_full::include_muons
bool include_muons
If true, include muons (default false)
Definition: eos_nse_full.h:185
o2scl::eos_nse_full::get_eos
o2scl::eos_had_temp_base & get_eos()
Get homogeneous matter EOS.
Definition: eos_nse_full.h:401
o2scl::eos_nse_full::solve_mun
virtual double solve_mun(double mun, dense_matter &dm)
Fix the baryon density by varying the neutron chemical potential.
o2scl::eos_nse_full::verbose
int verbose
Verbose parameter.
Definition: eos_nse_full.h:172
o2scl::eos_nse_full::ehtp
o2scl::eos_had_temp_base * ehtp
Nucleonic EOS (0 by default)
Definition: eos_nse_full.h:158

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