eos_nse_full.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2014-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_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
112 
113  /// Relativistic fermions with derivatives
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 Desc
203  */
204  virtual int bracket_mu_solve(double &mun_low, double &mun_high,
205  double &mup_low, double &mup_high,
206  dense_matter &dm);
207 
208  /** \brief Fix electron fraction by varying proton chemical
209  potential
210 
211  At some fixed values of <tt>dm.Ye</tt> and <tt>dm.nB</tt>,
212  given a value of \f$ \mu_p \f$, and given an initial bracket
213  for \f$ \mu_n \f$ (stored in <tt>mun_low</tt> and
214  <tt>mun_high</tt>), this function attempts to find the value
215  of \f$ \mu_n \f$ which ensures that the baryon density in
216  nuclei matches that in <tt>dm.nB</tt> using a bracketing
217  solver. It then returns the difference between the value of
218  the proton fraction in nuclei and the value in <tt>dm.Ye</tt>.
219 
220  If <tt>mun_low</tt> and <tt>mun_high</tt> do not bracket the
221  correct value of \f$ \mu_n \f$, this function attempts to
222  modify them to give a proper bracket for the root. The
223  finaly value of \f$ \mu_n \f$ is stored in <tt>dm.n.mu</tt>.
224 
225  Currently, the values of <tt>dm.n.n</tt> and <tt>dm.p.n</tt>
226  are ignored and set to zero.
227  */
228  double mup_for_Ye(double mup, double &mun_low,
229  double &mun_high, dense_matter &dm);
230 
231  /** \brief Fix the baryon density by varying the neutron
232  chemical potential
233 
234  Given a value of \f$ \mu_n \f$ (the value in <tt>dm.n.mu</tt>
235  is ignored), this function computes the baryon density
236  in nuclei and returns the difference between this value
237  and that stored in <tt>dm.nB</tt>.
238 
239  Currently, the values of <tt>dm.n.n</tt> and <tt>dm.p.n</tt>
240  are ignored and set to zero.
241  */
242  virtual double solve_mun(double mun, dense_matter &dm);
243 
244  /** \brief Compute the properties of matter from the densities,
245  not presuming equilibrium
246 
247  The values of <tt>dm.nB</tt> and <tt>dm.Ye</tt> are ignored
248  and unchanged by this function. The electron and muon density
249  are determined by charged neutrality and assuming their
250  chemical potentials are equal. Photons are always included.
251 
252  If the nuclear densities are all zero, then this just
253  returns nuclear matter with leptons and photons as
254  determined by charge neutrality.
255 
256  This function is designed to return non-zero values for
257  invalid configurations and can return the value
258  \ref invalid_config without calling the error handler,
259  independent of the value of \ref err_nonconv .
260 
261  Possible invalid configurations are:
262  - negative nucleon or nucleus densities, or
263  - proton radii larger than WS cell radii, i.e.
264  \f$ (0.08 - n_p) / (n_e+n_{\mu}-n_p) < 1 \f$ or
265  \f$ n_p > 0.08 \f$ .
266  */
267  virtual int calc_density_noneq(dense_matter &dm);
268 
269  /** \brief Compute the properties of matter from
270  neutron and proton densities, using the Saha equation
271 
272  If the parameter <tt>from_densities</tt> is true, then this
273  computes nucleonic matter using the neutron and proton
274  densities stored in <tt>dm.n.n</tt> and <tt>dm.p.n</tt>.
275  Otherwise, nucleonic matter is computed using the chemical
276  potential stored in <tt>dm.n.mu</tt> and <tt>dm.p.mu</tt>.
277  Either way, electrons are computed assuming their density is
278  given from \ref o2scl::dense_matter::nB and \ref
279  o2scl::dense_matter::Ye. Muons are added assuming their
280  chemical potential is equal to the electron chemical
281  potential. Finally, the Saha equation is used to determine the
282  nuclear chemical potentials and this gives the nuclear
283  densities.
284 
285  This function only works when \ref inc_prot_coul is
286  <tt>false</tt>.
287 
288  The values in \ref o2scl::dense_matter::nB and \ref
289  o2scl::dense_matter::Ye are unchanged by this function. Note
290  that, after this function completes, the value returned by
291  \ref o2scl::dense_matter::baryon_density() will not
292  necessarily be the same as that stored in \ref
293  o2scl::dense_matter::nB (and similarly for the electron
294  fraction).
295 
296  This function is designed to return non-zero values for
297  invalid configurations and can return the value
298  \ref invalid_config without calling the error handler,
299  independent of the value of \ref err_nonconv .
300 
301  Possible invalid configurations are:
302  - negative nucleon densities, or
303  - proton radii larger than WS cell radii, i.e.
304  \f$ (0.08 - n_p) / (n_e+n_{\mu}-n_p) < 1 \f$ or
305  \f$ n_p > 0.08 \f$ .
306  */
307  virtual int calc_density_fixnp(dense_matter &dm, bool from_densities=true);
308 
309  /** \brief Compute the free energy for a fixed composition
310  by minimization
311 
312  Given a fixed baryon density (dm.nB), electron fraction
313  (dm.Ye), temperature (dm.T), this minimizes the free energy
314  over the densities of the nuclei currently present in the
315  distribution. The neutron and proton drip densities are
316  determined by ensuring that the baryon density and electron
317  fraction are correctly reproduced. The function which
318  is minimized is \ref free_energy() .
319 
320  \note This function currently only performs a very simple
321  minimization and currently works in only limited
322  circumstances.
323  */
324  virtual int calc_density_by_min(dense_matter &dm);
325 
326  /** \brief Compute properties of matter for baryon density and
327  electron fraction using the Saha equation
328 
329  This function solves the function specified by \ref
330  solve_fixnp() using the current values of <tt>dm.n.n</tt> and
331  <tt>dm.p.n</tt> as initial guesses.
332  */
333  virtual int calc_density_saha(dense_matter &dm);
334 
335  /** \brief Output properties of a \ref o2scl::dense_matter object to
336  std::cout
337 
338  This function was particularly designed for comparing results
339  with \ref o2scl::eos_sn_base derived classes.
340 
341  If output level is 0, then just the basic quantities are
342  output without any information about the distribution. If
343  output_level is 1, then only about 10 nuclei in the
344  distribution are output, and if output_level is 2,
345  then all nuclei in the distribution are output.
346  */
347  virtual void output(dense_matter &dm, int output_level);
348 
349  /** \brief Adjust the particle densities to match specified
350  density and electron fraction
351 
352  This function attempts to match the nuclear and nucleon
353  densities so that the baryon density and electron fraction are
354  equal to those specified in \ref o2scl::dense_matter::nB and
355  \ref o2scl::dense_matter::Ye .
356  */
357  virtual int density_match(dense_matter &dm);
358 
359  /** \brief Relativistic fermions
360 
361  \comment
362  Must currently be public for tcan/ecn.cpp.
363  \endcomment
364  */
366 
367  /// \name Nuclei and nuclear masses
368  //@{
369  /// Compute nuclei in dense matter
371 
372  /** \brief Set nuclear mass formula
373  */
375  massp=&m;
376  return;
377  }
378 
379  /** \brief Set distribution of nuclei
380  */
381  void set_dist(std::vector<o2scl::nucleus> &dist) {
382  ad=&dist;
383  return;
384  }
385  //@}
386 
387  /// \name Nucleonic matter EOS
388  //@{
389  /** \brief Set homogeneous matter EOS
390  */
392  ehtp=&e;
393  return;
394  }
395 
396  /** \brief Get homogeneous matter EOS
397 
398  This function calls the error handler if no EOS has been set
399  */
401  if (ehtp==0) {
402  O2SCL_ERR2("Homogeneous matter EOS not specified in ",
403  "eos_nse_full::get_eos().",exc_efailed);
404  }
405  return *ehtp;
406  }
407 
408  /** \brief Return true if an EOS was specified
409  */
410  bool is_eos_set() {
411  if (ehtp==0) return false;
412  return true;
413  }
414  //@}
415 
416  /// \name Numerical methods
417  //@{
418  /// The default minimizer
420 
421  /// Default solver
423 
424  /// Lepton solver
426  //@}
427 
428 #ifdef O2SCL_NEVER_DEFINED
429 
430  /** \brief Desc
431 
432  This was intended to be a version of calc_density_by_min()
433  which optimized the composition, but it doesn't really work
434  yet.
435  */
436  int calc_density(dense_matter &dm);
437 
438 #endif
439 
440  };
441 
442 #ifndef DOXYGEN_NO_O2NS
443 }
444 #endif
445 
446 #endif
o2scl::fermion_deriv_rel snf
Relativistic fermions with derivatives.
Definition: eos_nse_full.h:114
std::vector< o2scl::nucleus > * ad
The full distribution of all nuclei to consider.
Definition: eos_nse_full.h:125
bool inc_prot_coul
If true, include dripped protons and neutrons in the nuclear mass (default true)
Definition: eos_nse_full.h:182
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...
void set_mass(o2scl::nucmass_densmat &m)
Set nuclear mass formula.
Definition: eos_nse_full.h:374
virtual double solve_mun(double mun, dense_matter &dm)
Fix the baryon density by varying the neutron chemical potential.
virtual double free_energy(const ubvector &n_nuc, dense_matter &dm)
Compute the free energy from a vector of densities of the nuclei.
o2scl::classical cla
Compute particle properties assuming classical thermodynamics.
Definition: eos_nse_full.h:111
int verbose
Verbose parameter.
Definition: eos_nse_full.h:172
virtual int calc_density_saha(dense_matter &dm)
Compute properties of matter for baryon density and electron fraction using the Saha equation...
void set_eos(o2scl::eos_had_temp_base &e)
Set homogeneous matter EOS.
Definition: eos_nse_full.h:391
mroot_hybrids def_mroot
Default solver.
Definition: eos_nse_full.h:422
root_cern def_root
Lepton solver.
Definition: eos_nse_full.h:425
A finite temperature hadronic EOS [abstract base].
Definition: eos_had_base.h:962
void set_dist(std::vector< o2scl::nucleus > &dist)
Set distribution of nuclei.
Definition: eos_nse_full.h:381
bool err_nonconv
If true, call the error handler if calc_density() does not converge (default true) ...
Definition: eos_nse_full.h:177
exc_efailed
virtual int calc_density_by_min(dense_matter &dm)
Compute the free energy for a fixed composition by minimization.
o2scl::nucmass_densmat * massp
Mass formula (points to nuc_dens by default)
Definition: eos_nse_full.h:117
EOS for nuclear statistical equilibrium with interactions.
Definition: eos_nse_full.h:84
o2scl::eos_had_temp_base * ehtp
Nucleonic EOS (0 by default)
Definition: eos_nse_full.h:158
virtual int density_match(dense_matter &dm)
Adjust the particle densities to match specified density and electron fraction.
o2scl::nucmass_densmat nuc_dens
Compute nuclei in dense matter.
Definition: eos_nse_full.h:370
o2scl::fermion_rel relf
Relativistic fermions.
Definition: eos_nse_full.h:365
static const int invalid_config
The integer return value which indicates an invalid configuration.
Definition: eos_nse_full.h:167
double mup_for_Ye(double mup, double &mun_low, double &mun_high, dense_matter &dm)
Fix electron fraction by varying proton chemical potential.
#define O2SCL_ERR2(d, d2, n)
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_had_temp_base & get_eos()
Get homogeneous matter EOS.
Definition: eos_nse_full.h:400
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::mmin_simp2 def_mmin
The default minimizer.
Definition: eos_nse_full.h:419
bool include_muons
If true, include muons (default false)
Definition: eos_nse_full.h:185
virtual void output(dense_matter &dm, int output_level)
Output properties of a o2scl::dense_matter object to std::cout.
virtual int bracket_mu_solve(double &mun_low, double &mun_high, double &mup_low, double &mup_high, dense_matter &dm)
Desc.
bool inc_lept_phot
If true, include electrons and photons (default true)
Definition: eos_nse_full.h:108
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...
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...
bool is_eos_set()
Return true if an EOS was specified.
Definition: eos_nse_full.h:410
virtual int calc_density_noneq(dense_matter &dm)
Compute the properties of matter from the densities, not presuming equilibrium.

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