part.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-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 #ifndef O2SCL_PART_H
24 #define O2SCL_PART_H
25 
26 #include <string>
27 #include <iostream>
28 #include <cmath>
29 #include <o2scl/constants.h>
30 #include <o2scl/inte.h>
31 #include <o2scl/funct.h>
32 #include <o2scl/mroot.h>
33 
34 // To get directories for calibrate function
35 #include <o2scl/lib_settings.h>
36 // To read tables in calibrate function
37 #include <o2scl/table.h>
38 #include <o2scl/hdf_file.h>
39 #include <o2scl/hdf_io.h>
40 
41 /** \file part.h
42  \brief File defining \ref o2scl::thermo_tl and \ref o2scl::part_tl
43 */
44 
45 #ifndef DOXYGEN_NO_O2NS
46 namespace o2scl {
47 #endif
48 
49  /** \brief A class containing three thermodynamical variables (energy
50  density, pressure, entropy density)
51  */
52  template<class fp_t=double> class thermo_tl {
53 
54  public:
55 
56  /// Pressure
57  fp_t pr;
58  /// Energy density
59  fp_t ed;
60  /// Entropy density
61  fp_t en;
62 
63  /// Return string denoting type ("thermo")
64  const char *type() { return "thermo"; }
65 
66  // Default constructor
67  thermo_tl() {
68  }
69 
70  /// Copy constructor
71  thermo_tl(const thermo_tl &t) {
72  ed=t.ed;
73  pr=t.pr;
74  en=t.en;
75  }
76 
77  /// Copy construction with operator=()
79  if (this!=&t) {
80  ed=t.ed;
81  pr=t.pr;
82  en=t.en;
83  }
84  return *this;
85  }
86 
87  };
88 
89  /** \brief Double-precision thermodynamics object
90  */
92 
93  /** \brief Addition operator
94  */
95  extern thermo operator+(const thermo &left, const thermo &right);
96 
97  /** \brief Subtraction operator
98  */
99  extern thermo operator-(const thermo &left, const thermo &right);
100 
101  /** \brief Particle base class
102  */
103  template<class fp_t=double> class part_tl {
104 
105  public:
106 
107  /// Degeneracy (e.g. spin and color if applicable)
108  fp_t g;
109  /// Mass
110  fp_t m;
111  /// Number density
112  fp_t n;
113  /// Energy density
114  fp_t ed;
115  /// Pressure
116  fp_t pr;
117  /// Chemical potential
118  fp_t mu;
119  /// Entropy density
120  fp_t en;
121  /// Effective mass (Dirac unless otherwise specified)
122  fp_t ms;
123  /// Effective chemical potential
124  fp_t nu;
125  /** \brief If true, include the mass in the energy
126  density and chemical potential (default true)
127  */
129  /// True if the particle is non-interacting (default true)
131 
132  /// Copy constructor
133  part_tl(const part_tl &p) {
134  g=p.g;
135  m=p.m;
136  ms=p.ms;
137  n=p.n;
138  ed=p.ed;
139  pr=p.pr;
140  mu=p.mu;
141  en=p.en;
142  nu=p.nu;
145  }
146 
147  /// Copy construction with operator=()
148  part_tl &operator=(const part_tl &p) {
149  if (this!=&p) {
150  g=p.g;
151  m=p.m;
152  ms=p.ms;
153  n=p.n;
154  ed=p.ed;
155  pr=p.pr;
156  mu=p.mu;
157  en=p.en;
158  nu=p.nu;
161  }
162  return *this;
163  }
164 
165  /// Make a particle of mass \c mass and degeneracy \c dof.
166  part_tl(fp_t mass=0.0, fp_t dof=0.0) {
167  m=mass;
168  ms=mass;
169  g=dof;
170 
171  non_interacting=true;
172  inc_rest_mass=true;
173  }
174 
175  virtual ~part_tl() {
176  }
177 
178  /// Set the mass \c mass and degeneracy \c dof.
179  virtual void init(fp_t mass, fp_t dof) {
180  m=mass;
181  ms=mass;
182  g=dof;
183  return;
184  }
185 
186  /** \brief Make \c ap an anti-particle with the same mass
187  and degeneracy
188 
189  This sets the \ref m, \ref g, \ref ms, \ref inc_rest_mass
190  and \ref non_interacting fields of \c ap equal to that
191  of the current object. If \ref inc_rest_mass is true,
192  then it sets
193  \f[
194  \mu_{\mathrm{anti}} = - \mu
195  \qquad\mathrm{and}\qquad
196  \nu_{\mathrm{anti}} = - \nu
197  \f]
198  and if \ref inc_rest_mass is false, it sets
199  \f[
200  \mu_{\mathrm{anti}} = - \mu - 2 m
201  \qquad\mathrm{and}\qquad
202  \nu_{\mathrm{anti}} = - \nu - 2 m
203  \f]
204  */
205  virtual void anti(part_tl &ax) {
206  ax.g=g;
207  ax.m=m;
208  ax.ms=ms;
211  if (inc_rest_mass) {
212  ax.nu=-nu;
213  ax.mu=-mu;
214  } else {
215  ax.nu=-nu-2.0*m;
216  ax.mu=-mu-2.0*m;
217  }
218  return;
219  }
220 
221  /// Return string denoting type ("part_tl")
222  virtual const char *type() { return "part_tl"; }
223 
224  };
225 
226  /** \brief Double-precision thermodynamics object
227  */
229 
230  /** \brief Addition operator
231  */
232  extern thermo operator+(const thermo &left, const part &right);
233 
234  /** \brief Subtraction operator
235  */
236  extern thermo operator-(const thermo &left, const part &right);
237 
238  /** \brief Object to organize calibration of particle classes
239  to results stored in a table
240  */
242 
243  public:
244 
245  /** \brief Set mass and flags from mot, T, and the index k
246  */
247  template<class part_t>
248  void set_mass_flags(part_t &p, double mot, double T, size_t k) {
249  if (k>=2) {
250  p.non_interacting=false;
251  p.ms=mot*T;
252  p.m=p.ms*1.5;
253  } else {
254  p.non_interacting=true;
255  p.m=mot*T;
256  p.ms=0.0;
257  }
258  if (k%2==0) {
259  p.inc_rest_mass=true;
260  } else {
261  p.inc_rest_mass=false;
262  }
263  return;
264  }
265 
266  /** \brief Set chemical potential from psi, T, the index k,
267  and the flag nr_mode
268  */
269  template<class part_t>
270  void set_chem_pot(part_t &p, double psi, double T, size_t k,
271  bool nr_mode) {
272  if (k%2==0) {
273  if (k>=2) {
274  if (nr_mode) {
275  p.nu=T*psi+p.m;
276  } else {
277  p.nu=T*psi+p.ms;
278  }
279  p.mu=0.0;
280  } else {
281  p.nu=0.0;
282  p.mu=T*psi+p.m;
283  }
284  } else {
285  if (k>=2) {
286  if (nr_mode) {
287  p.nu=T*psi;
288  } else {
289  p.nu=T*psi-p.m+p.ms;
290  }
291  p.mu=0.0;
292  } else {
293  p.nu=0.0;
294  p.mu=T*psi;
295  }
296  }
297  return;
298  }
299 
300  /** \brief Check the density against the exact result
301  and update
302  */
303  template<class part1_t, class part2_t, class part3_t>
304  void check_density(part1_t &p, part2_t &exact, part3_t &bad, size_t k,
305  double T, double mot, double psi,
306  double &mu_bad, double &m_bad,
307  double &T_bad, double &mot_bad, double &psi_bad,
308  double &ret_local) {
309  if (fabs((p.n-exact.n)/exact.n)>bad.n) {
310  bad.n=fabs((p.n-exact.n)/exact.n);
311  if (bad.n>ret_local) {
312  if (k>=2) {
313  mu_bad=p.nu;
314  m_bad=p.ms;
315  } else {
316  mu_bad=p.mu;
317  m_bad=p.m;
318  }
319  T_bad=T;
320  mot_bad=mot;
321  psi_bad=psi;
322  ret_local=bad.n;
323  }
324  }
325  return;
326  }
327 
328  /** \brief Check the chemical potential against the exact result
329  and update
330  */
331  template<class part1_t, class part2_t, class part3_t>
332  void check_chem_pot(part1_t &p, part2_t &exact, part3_t &bad, size_t k,
333  double T, double mot, double psi,
334  double &mu_bad, double &m_bad,
335  double &T_bad, double &mot_bad, double &psi_bad,
336  double &ret_local) {
337  if (k>=2) {
338  if (fabs((p.nu-exact.nu)/exact.nu)>bad.mu) {
339  bad.mu=fabs((p.nu-exact.nu)/exact.nu);
340  if (bad.mu>ret_local) {
341  mu_bad=p.nu;
342  m_bad=p.ms;
343  T_bad=T;
344  mot_bad=mot;
345  psi_bad=psi;
346  ret_local=bad.n;
347  }
348  }
349  } else {
350  if (fabs((p.mu-exact.mu)/exact.mu)>bad.mu) {
351  bad.mu=fabs((p.mu-exact.mu)/exact.mu);
352  if (bad.mu>ret_local) {
353  mu_bad=p.mu;
354  m_bad=p.m;
355  T_bad=T;
356  mot_bad=mot;
357  psi_bad=psi;
358  ret_local=bad.n;
359  }
360  }
361  }
362  return;
363  }
364 
365  /** \brief Check the energy density, pressure, and entropy against
366  the exact result and update
367  */
368  template<class part1_t, class part2_t, class part3_t>
369  void check_eps(part1_t &p, part2_t &exact, part3_t &bad, size_t k,
370  double T, double mot, double psi,
371  double &mu_bad, double &m_bad,
372  double &T_bad, double &mot_bad, double &psi_bad,
373  double &ret_local) {
374  if (fabs((p.ed-exact.ed)/exact.ed)>bad.ed) {
375  bad.ed=fabs((p.ed-exact.ed)/exact.ed);
376  if (bad.ed>ret_local) {
377  if (k>=2) {
378  mu_bad=p.nu;
379  m_bad=p.ms;
380  } else {
381  mu_bad=p.mu;
382  m_bad=p.m;
383  }
384  T_bad=T;
385  mot_bad=mot;
386  psi_bad=psi;
387  ret_local=bad.ed;
388  }
389  }
390  if (fabs((p.pr-exact.pr)/exact.pr)>bad.pr) {
391  bad.pr=fabs((p.pr-exact.pr)/exact.pr);
392  if (bad.pr>ret_local) {
393  if (k>=2) {
394  mu_bad=p.nu;
395  m_bad=p.ms;
396  } else {
397  mu_bad=p.mu;
398  m_bad=p.m;
399  }
400  T_bad=T;
401  mot_bad=mot;
402  psi_bad=psi;
403  ret_local=bad.pr;
404  }
405  }
406  if (fabs((p.en-exact.en)/exact.en)>bad.en) {
407  bad.en=fabs((p.en-exact.en)/exact.en);
408  if (bad.en>ret_local) {
409  if (k>=2) {
410  mu_bad=p.nu;
411  m_bad=p.ms;
412  } else {
413  mu_bad=p.mu;
414  m_bad=p.m;
415  }
416  T_bad=T;
417  mot_bad=mot;
418  psi_bad=psi;
419  ret_local=bad.en;
420  }
421  }
422  return;
423  }
424 
425  /** \brief Check the three derivatives against
426  the exact result and update
427  */
428  template<class part1_t, class part2_t, class part3_t>
429  void check_derivs(part1_t &p, part2_t &exact, part3_t &bad, size_t k,
430  double T, double mot, double psi,
431  double &mu_bad, double &m_bad,
432  double &T_bad, double &mot_bad, double &psi_bad,
433  double &ret_local) {
434  if (fabs((p.dndT-exact.dndT)/exact.dndT)>bad.dndT) {
435  bad.dndT=fabs((p.dndT-exact.dndT)/exact.dndT);
436  if (bad.dndT>ret_local) {
437  if (k>=2) {
438  mu_bad=p.nu;
439  m_bad=p.ms;
440  } else {
441  mu_bad=p.mu;
442  m_bad=p.m;
443  }
444  T_bad=T;
445  mot_bad=mot;
446  psi_bad=psi;
447  ret_local=bad.dndT;
448  }
449  }
450  if (fabs((p.dndmu-exact.dndmu)/exact.dndmu)>bad.dndmu) {
451  bad.dndmu=fabs((p.dndmu-exact.dndmu)/exact.dndmu);
452  if (bad.dndmu>ret_local) {
453  if (k>=2) {
454  mu_bad=p.nu;
455  m_bad=p.ms;
456  } else {
457  mu_bad=p.mu;
458  m_bad=p.m;
459  }
460  T_bad=T;
461  mot_bad=mot;
462  psi_bad=psi;
463  ret_local=bad.dndmu;
464  }
465  }
466  if (fabs((p.dsdT-exact.dsdT)/exact.dsdT)>bad.dsdT) {
467  bad.dsdT=fabs((p.dsdT-exact.dsdT)/exact.dsdT);
468  if (bad.dsdT>ret_local) {
469  if (k>=2) {
470  mu_bad=p.nu;
471  m_bad=p.ms;
472  } else {
473  mu_bad=p.mu;
474  m_bad=p.m;
475  }
476  T_bad=T;
477  mot_bad=mot;
478  psi_bad=psi;
479  ret_local=bad.dsdT;
480  }
481  }
482  return;
483  }
484 
485  /** \brief Calibrate a particle thermodynamics class with results
486  stored in a table
487 
488  This compares the approximation to the exact results using
489  calc_density(), calc_mu(), pair_density() and pair_mu(). It
490  tries each function twelve times. It tries three different
491  temperatures, setting both <tt>inc_rest_mass</tt> and
492  <tt>non_interacting</tt> equal to <tt>true</tt> and
493  <tt>false</tt>.
494 
495  The <tt>verbose</tt> parameter controls the amount of output.
496 
497  \future Also calibrate massless fermions?
498  */
499  template<class part_t, class thermo_t>
500  double part_calibrate(part_t &p, thermo_t &th, bool test_pair,
501  std::string file, bool nr_mode=false,
502  int verbose=0, bool external=false) {
503 
504  double ret=0;
505 
506  // ----------------------------------------------------------------
507  // Will return to these original values afterwards
508 
509  part_t orig=p;
510 
511  // ----------------------------------------------------------------
512  // Read data file
513 
514  std::string fname;
515  if (external==false) {
516  fname=o2scl_settings.get_data_dir()+file;
517  } else {
518  fname=file;
519  }
520 
521  if (verbose>1) {
522  std::cout << "In part_calibrate(), loading file named\n\t'"
523  << fname << "'.\n" << std::endl;
524  }
525  o2scl::table<> tab;
527  hf.open(fname);
528  std::string name;
529 #ifndef O2SCL_NO_HDF_INPUT
530  hdf_input(hf,tab,name);
531 #endif
532  hf.close();
533 
534  if (tab.get_nlines()==0) {
535  std::string str="Failed to load data from file '"+fname+
536  "' in part_calibrate(). Bad filename?";
537  O2SCL_ERR(str.c_str(),exc_efilenotfound);
538  }
539 
540  if (!tab.is_column("ed")) {
541  tab.function_column("ed_mot*mot","ed");
542  if (test_pair) {
543  tab.function_column("pair_ed_mot*mot","pair_ed");
544  }
545  }
546 
547  // ----------------------------------------------------------------
548 
549  p.g=2.0;
550 
551  size_t cnt=0;
552  part bad, dev, exact;
553  double m_bad=0.0, mu_bad=0.0, T_bad=0.0, mot_bad=0.0, psi_bad=0.0;
554  p.non_interacting=true;
555 
556  // ----------------------------------------------------------------
557  // First pass, test calc_mu()
558 
559  // k=0,2 are with rest mass, k=1,3 are without
560  // k=0,1 are non-interacting, k=2,3 are interacting
561  for(size_t k=0;k<4;k++) {
562 
563  double ret_local=0.0;
564 
565  // Initialize storage
566  dev.n=0.0; dev.ed=0.0; dev.pr=0.0; dev.en=0.0;
567  bad.n=0.0; bad.ed=0.0; bad.pr=0.0; bad.en=0.0;
568 
569  // Temperature loop
570  for(double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
571 
572  // Loop over each point in the data file
573  for(size_t i=0;i<tab.get_nlines();i++) {
574 
575  double mot=tab.get("mot",i);
576  double psi=tab.get("psi",i);
577  exact.n=tab.get("n",i);
578  exact.ed=tab.get("ed",i);
579  exact.pr=tab.get("pr",i);
580  exact.en=tab.get("en",i);
581 
582  set_mass_flags(p,mot,T,k);
583  set_chem_pot(p,psi,T,k,nr_mode);
584 
585  th.calc_mu(p,T);
586 
587  exact.n*=pow(T,3.0);
588  if (nr_mode) {
589  if (k%2==0) {
590  exact.ed=exact.ed*pow(T,4.0)+exact.n*p.m;
591  } else {
592  exact.ed=exact.ed*pow(T,4.0);
593  }
594  } else {
595  if (k%2==0) {
596  exact.ed*=pow(T,4.0);
597  } else {
598  exact.ed=exact.ed*pow(T,4.0)-exact.n*p.m;
599  }
600  }
601  exact.pr*=pow(T,4.0);
602  exact.en*=pow(T,3.0);
603 
604  dev.n+=fabs((p.n-exact.n)/exact.n);
605  dev.ed+=fabs((p.ed-exact.ed)/exact.ed);
606  dev.pr+=fabs((p.pr-exact.pr)/exact.pr);
607  dev.en+=fabs((p.en-exact.en)/exact.en);
608 
609  cnt++;
610 
611  check_density<part_t>(p,exact,bad,k,T,mot,psi,mu_bad,m_bad,T_bad,
612  mot_bad,psi_bad,ret_local);
613 
614  check_eps<part_t>(p,exact,bad,k,T,mot,psi,mu_bad,m_bad,T_bad,
615  mot_bad,psi_bad,ret_local);
616 
617  if (verbose>1) {
618  std::cout.precision(5);
619  if (k>=2) {
620  std::cout << "T,ms,nu,psi,mot: " << T << " "
621  << p.ms << " " << p.nu << " "
622  << psi << " " << mot << std::endl;
623  } else {
624  std::cout << "T,m,mu,psi,mot: " << T << " "
625  << p.m << " " << p.mu << " "
626  << psi << " " << mot << std::endl;
627  }
628  std::cout.precision(5);
629  std::cout << "n,ed,pr,en: " << std::endl;
630  std::cout << "approx: " << p.n << " " << p.ed << " "
631  << p.pr << " " << p.en << std::endl;
632  std::cout << "exact : " << exact.n << " " << exact.ed << " "
633  << exact.pr << " " << exact.en << std::endl;
634  std::cout << "bad : " << bad.n << " " << bad.ed << " "
635  << bad.pr << " " << bad.en << std::endl;
636  std::cout << "ret_local,ret: " << ret_local << " "
637  << ret << std::endl;
638  std::cout << std::endl;
639  if (verbose>2) {
640  char ch;
641  std::cin >> ch;
642  }
643  }
644 
645  if (ret_local>ret) {
646  ret=ret_local;
647  }
648 
649  // End of loop over points in data file
650  }
651  // End of temperature loop
652  }
653 
654  dev.n/=cnt;
655  dev.ed/=cnt;
656  dev.pr/=cnt;
657  dev.en/=cnt;
658 
659  if (verbose>0) {
660  if (k==0) {
661  std::cout << "Function calc_mu(), include rest mass"
662  << std::endl;
663  } else if (k==1) {
664  std::cout << "Function calc_mu(), without rest mass"
665  << std::endl;
666  } else if (k==2) {
667  std::cout << "Function calc_mu(), include rest mass, "
668  << "interacting" << std::endl;
669  } else {
670  std::cout << "Function calc_mu(), without rest mass, "
671  << "interacting" << std::endl;
672  }
673 
674  std::cout << "Average performance: " << std::endl;
675  std::cout << "n: " << dev.n << " ed: " << dev.ed << " pr: "
676  << dev.pr << " en: " << dev.en << std::endl;
677  std::cout << "Worst case: " << std::endl;
678  std::cout << "n: " << bad.n << " ed: " << bad.ed << " pr: "
679  << bad.pr << " en: " << bad.en << std::endl;
680  std::cout << "mu: " << mu_bad << " m: " << m_bad
681  << " T: " << T_bad << " mot: " << mot_bad
682  << "\n\tpsi: " << psi_bad << std::endl;
683  std::cout << "ret_local,ret: " << ret_local << " "
684  << ret << std::endl;
685  std::cout << std::endl;
686  if (verbose>2) {
687  char ch;
688  std::cin >> ch;
689  }
690  }
691 
692  // Reset p.non_interacting
693  p.non_interacting=true;
694 
695  // End of k loop
696  }
697 
698  // ----------------------------------------------------------------
699  // Second pass, test calc_density()
700 
701  // k=0,2 are with rest mass, k=1,3 are without
702  // k=0,1 are non-interacting, k=2,3 are interacting
703  for(size_t k=0;k<4;k++) {
704 
705  double ret_local=0.0;
706 
707  // Initialize storage
708  dev.mu=0.0; dev.ed=0.0; dev.pr=0.0; dev.en=0.0;
709  bad.mu=0.0; bad.ed=0.0; bad.pr=0.0; bad.en=0.0;
710 
711  // Temperature loop
712  for(double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
713 
714  // Loop over each point in the data file
715  for(size_t i=0;i<tab.get_nlines();i++) {
716 
717  double mot=tab.get("mot",i);
718  double psi=tab.get("psi",i);
719  p.n=tab.get("n",i)*pow(T,3.0);
720  exact.ed=tab.get("ed",i);
721  exact.pr=tab.get("pr",i);
722  exact.en=tab.get("en",i);
723 
724  set_mass_flags(p,mot,T,k);
725  set_chem_pot(exact,psi,T,k,nr_mode);
726 
727  exact.n=p.n;
728  if (nr_mode) {
729  if (k%2==0) {
730  exact.ed=exact.ed*pow(T,4.0)+exact.n*p.m;
731  } else {
732  exact.ed=exact.ed*pow(T,4.0);
733  }
734  } else {
735  if (k%2==0) {
736  exact.ed*=pow(T,4.0);
737  } else {
738  exact.ed=exact.ed*pow(T,4.0)-exact.n*p.m;
739  }
740  }
741  exact.pr*=pow(T,4.0);
742  exact.en*=pow(T,3.0);
743 
744  // Give it a guess for the chemical potential
745  if (k>=2) {
746  p.nu=p.m;
747  } else {
748  p.mu=p.m;
749  }
750 
751  th.calc_density(p,T);
752 
753  if (k>=2) {
754  dev.nu+=fabs((p.nu-exact.nu)/exact.nu);
755  } else {
756  dev.mu+=fabs((p.mu-exact.mu)/exact.mu);
757  }
758  dev.ed+=fabs((p.ed-exact.ed)/exact.ed);
759  dev.pr+=fabs((p.pr-exact.pr)/exact.pr);
760  dev.en+=fabs((p.en-exact.en)/exact.en);
761 
762  cnt++;
763 
764  check_chem_pot<part_t>(p,exact,bad,k,T,mot,psi,mu_bad,m_bad,T_bad,
765  mot_bad,psi_bad,ret_local);
766 
767  check_eps<part_t>(p,exact,bad,k,T,mot,psi,mu_bad,m_bad,T_bad,
768  mot_bad,psi_bad,ret_local);
769 
770  if (verbose>1) {
771  std::cout.precision(5);
772  if (k>=2) {
773  std::cout << "T,ms,n,psi,mot: " << T << " "
774  << p.ms << " " << p.n << " "
775  << psi << " " << mot << std::endl;
776  } else {
777  std::cout << "T,m,n,psi,mot: " << T << " "
778  << p.m << " " << p.n << " "
779  << psi << " " << mot << std::endl;
780  }
781  std::cout.precision(6);
782  if (k>=2) {
783  std::cout << "nu,ed,pr,en: " << std::endl;
784  std::cout << "approx: " << p.nu << " " << p.ed << " "
785  << p.pr << " " << p.en << std::endl;
786  std::cout << "exact : " << exact.nu << " " << exact.ed << " "
787  << exact.pr << " " << exact.en << std::endl;
788  } else {
789  std::cout << "mu,ed,pr,en: " << std::endl;
790  std::cout << "approx: " << p.mu << " " << p.ed << " "
791  << p.pr << " " << p.en << std::endl;
792  std::cout << "exact : " << exact.mu << " " << exact.ed << " "
793  << exact.pr << " " << exact.en << std::endl;
794  }
795  std::cout << "bad : " << bad.mu << " " << bad.ed << " "
796  << bad.pr << " " << bad.en << std::endl;
797  std::cout << "ret_local,ret: " << ret_local << " "
798  << ret << std::endl;
799  std::cout << std::endl;
800  if (verbose>2) {
801  char ch;
802  std::cin >> ch;
803  }
804  }
805 
806  if (ret_local>ret) {
807  ret=ret_local;
808  }
809 
810  // End of loop over points in data file
811  }
812  // End of temperature loop
813  }
814 
815  dev.mu/=cnt;
816  dev.ed/=cnt;
817  dev.pr/=cnt;
818  dev.en/=cnt;
819 
820  if (verbose>0) {
821  if (k==0) {
822  std::cout << "Function calc_density(), include rest mass"
823  << std::endl;
824  } else if (k==1) {
825  std::cout << "Function calc_density(), without rest mass"
826  << std::endl;
827  } else if (k==2) {
828  std::cout << "Function calc_density(), include rest mass, "
829  << "interacting" << std::endl;
830  } else {
831  std::cout << "Function calc_density(), without rest mass, "
832  << "interacting" << std::endl;
833  }
834 
835  std::cout << "Average performance: " << std::endl;
836  std::cout << "mu: " << dev.mu << " ed: " << dev.ed << " pr: "
837  << dev.pr << " en: " << dev.en << std::endl;
838  std::cout << "Worst case: " << std::endl;
839  std::cout << "mu: " << bad.mu << " ed: " << bad.ed << " pr: "
840  << bad.pr << " en: " << bad.en << std::endl;
841  std::cout << "mu: " << mu_bad << " m: " << m_bad
842  << " T: " << T_bad << " mot: " << mot_bad
843  << "\n\tpsi: " << psi_bad << std::endl;
844  std::cout << "ret_local,ret: " << ret_local << " "
845  << ret << std::endl;
846  std::cout << std::endl;
847  if (verbose>2) {
848  char ch;
849  std::cin >> ch;
850  }
851  }
852 
853  // End of k loop
854  }
855 
856  if (test_pair) {
857 
858  // ----------------------------------------------------------------
859  // Third pass, test pair_mu()
860 
861  // k=0,2 are with rest mass, k=1,3 are without
862  // k=0,1 are non-interacting, k=2,3 are interacting
863  for(size_t k=0;k<4;k++) {
864 
865  double ret_local=0.0;
866 
867  // Initialize storage
868  dev.n=0.0; dev.ed=0.0; dev.pr=0.0; dev.en=0.0;
869  bad.n=0.0; bad.ed=0.0; bad.pr=0.0; bad.en=0.0;
870 
871  // Temperature loop
872  for(double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
873 
874  // Loop over each point in the data file
875  for(size_t i=0;i<tab.get_nlines();i++) {
876 
877  double mot=tab.get("mot",i);
878  double psi=tab.get("psi",i);
879  exact.n=tab.get("pair_n",i);
880  exact.ed=tab.get("pair_ed",i);
881  exact.pr=tab.get("pair_pr",i);
882  exact.en=tab.get("pair_en",i);
883 
884  set_mass_flags(p,mot,T,k);
885  set_chem_pot(p,psi,T,k,nr_mode);
886 
887  th.pair_mu(p,T);
888 
889  exact.n*=pow(T,3.0);
890  if (k%2==0) {
891  exact.ed*=pow(T,4.0);
892  } else {
893  exact.ed=exact.ed*pow(T,4.0)-exact.n*p.m;
894  }
895  exact.pr*=pow(T,4.0);
896  exact.en*=pow(T,3.0);
897 
898  dev.n+=fabs((p.n-exact.n)/exact.n);
899  dev.ed+=fabs((p.ed-exact.ed)/exact.ed);
900  dev.pr+=fabs((p.pr-exact.pr)/exact.pr);
901  dev.en+=fabs((p.en-exact.en)/exact.en);
902 
903  cnt++;
904 
905  check_density<part_t>(p,exact,bad,k,T,mot,psi,mu_bad,m_bad,T_bad,
906  mot_bad,psi_bad,ret_local);
907 
908  check_eps<part_t>(p,exact,bad,k,T,mot,psi,mu_bad,m_bad,T_bad,
909  mot_bad,psi_bad,ret_local);
910 
911  if (verbose>1) {
912  std::cout.precision(5);
913  std::cout << "T,m,mu,psi,mot: " << T << " " << p.m << " "
914  << p.mu << " " << psi << " " << mot << std::endl;
915  std::cout.precision(6);
916  std::cout << i << " " << exact.m << " " << exact.ms << " "
917  << exact.mu << " " << exact.nu << std::endl;
918  std::cout << "n,ed,pr,en: " << std::endl;
919  std::cout << "approx: " << p.n << " " << p.ed << " "
920  << p.pr << " " << p.en << std::endl;
921  std::cout << "exact : " << exact.n << " " << exact.ed << " "
922  << exact.pr << " " << exact.en << std::endl;
923  std::cout << "bad : " << bad.n << " " << bad.ed << " "
924  << bad.pr << " " << bad.en << std::endl;
925  std::cout << "ret_local,ret: " << ret_local << " "
926  << ret << std::endl;
927  std::cout << std::endl;
928  if (verbose>2) {
929  char ch;
930  std::cin >> ch;
931  }
932  }
933 
934  if (ret_local>ret) {
935  ret=ret_local;
936  }
937 
938  // End of loop over points in data file
939  }
940  // End of temperature loop
941  }
942 
943  dev.n/=cnt;
944  dev.ed/=cnt;
945  dev.pr/=cnt;
946  dev.en/=cnt;
947 
948  if (verbose>0) {
949  if (k==0) {
950  std::cout << "Function pair_mu(), include rest mass"
951  << std::endl;
952  } else if (k==1) {
953  std::cout << "Function pair_mu(), without rest mass"
954  << std::endl;
955  } else if (k==2) {
956  std::cout << "Function pair_mu(), include rest mass, "
957  << "interacting" << std::endl;
958  } else {
959  std::cout << "Function pair_mu(), without rest mass, "
960  << "interacting" << std::endl;
961  }
962 
963  std::cout << "Average performance: " << std::endl;
964  std::cout << "n: " << dev.n << " ed: " << dev.ed << " pr: "
965  << dev.pr << " en: " << dev.en << std::endl;
966  std::cout << "Worst case: " << std::endl;
967  std::cout << "n: " << bad.n << " ed: " << bad.ed << " pr: "
968  << bad.pr << " en: " << bad.en << std::endl;
969  std::cout << "mu: " << mu_bad << " m: " << m_bad
970  << " T: " << T_bad << " mot: " << mot_bad
971  << "\n\tpsi: " << psi_bad << std::endl;
972  std::cout << "ret_local,ret: " << ret_local << " "
973  << ret << std::endl;
974  std::cout << std::endl;
975  if (verbose>2) {
976  char ch;
977  std::cin >> ch;
978  }
979  }
980 
981  // End of k loop
982  }
983 
984  // ----------------------------------------------------------------
985  // Fourth pass, test pair_density()
986 
987  // k=0,2 are with rest mass, k=1,3 are without
988  // k=0,1 are non-interacting, k=2,3 are interacting
989  for(size_t k=0;k<4;k++) {
990 
991  double ret_local=0.0;
992 
993  // Initialize storage
994  dev.mu=0.0; dev.ed=0.0; dev.pr=0.0; dev.en=0.0;
995  bad.mu=0.0; bad.ed=0.0; bad.pr=0.0; bad.en=0.0;
996 
997  // Temperature loop
998  for(double T=1.0e-2;T<=1.001e2;T*=1.0e2) {
999 
1000  // Loop over each point in the data file
1001  for(size_t i=0;i<tab.get_nlines();i++) {
1002 
1003  double mot=tab.get("mot",i);
1004  double psi=tab.get("psi",i);
1005  p.n=tab.get("pair_n",i)*pow(T,3.0);
1006  exact.ed=tab.get("pair_ed",i);
1007  exact.pr=tab.get("pair_pr",i);
1008  exact.en=tab.get("pair_en",i);
1009 
1010  set_mass_flags(p,mot,T,k);
1011  set_chem_pot(exact,psi,T,k,nr_mode);
1012 
1013  exact.n=p.n;
1014  if (k%2==0) {
1015  exact.ed*=pow(T,4.0);
1016  } else {
1017  exact.ed=exact.ed*pow(T,4.0)-exact.n*p.m;
1018  }
1019  exact.pr*=pow(T,4.0);
1020  exact.en*=pow(T,3.0);
1021 
1022  // Give it a guess for the chemical potential
1023  p.mu=p.m;
1024 
1025  th.pair_density(p,T);
1026 
1027  if (k>=2) {
1028  dev.nu+=fabs((p.nu-exact.nu)/exact.nu);
1029  } else {
1030  dev.mu+=fabs((p.mu-exact.mu)/exact.mu);
1031  }
1032  dev.ed+=fabs((p.ed-exact.ed)/exact.ed);
1033  dev.pr+=fabs((p.pr-exact.pr)/exact.pr);
1034  dev.en+=fabs((p.en-exact.en)/exact.en);
1035 
1036  cnt++;
1037 
1038  check_chem_pot<part_t>(p,exact,bad,k,T,mot,psi,mu_bad,m_bad,T_bad,
1039  mot_bad,psi_bad,ret_local);
1040 
1041  check_eps<part_t>(p,exact,bad,k,T,mot,psi,mu_bad,m_bad,T_bad,
1042  mot_bad,psi_bad,ret_local);
1043 
1044 
1045  if (verbose>1) {
1046  std::cout.precision(5);
1047  if (k>=2) {
1048  std::cout << "T,ms,n,psi,mot: " << T << " " << p.ms << " "
1049  << p.n << " " << psi << " " << mot << std::endl;
1050  } else {
1051  std::cout << "T,m,n,psi,mot: " << T << " " << p.m << " "
1052  << p.n << " " << psi << " " << mot << std::endl;
1053  }
1054  std::cout.precision(6);
1055  if (k>=2) {
1056  std::cout << "nu,ed,pr,en: " << std::endl;
1057  std::cout << "approx: " << p.nu << " " << p.ed << " "
1058  << p.pr << " " << p.en << std::endl;
1059  std::cout << "exact : " << exact.nu << " " << exact.ed << " "
1060  << exact.pr << " " << exact.en << std::endl;
1061  } else {
1062  std::cout << "mu,ed,pr,en: " << std::endl;
1063  std::cout << "approx: " << p.mu << " " << p.ed << " "
1064  << p.pr << " " << p.en << std::endl;
1065  std::cout << "exact : " << exact.mu << " " << exact.ed << " "
1066  << exact.pr << " " << exact.en << std::endl;
1067  }
1068  std::cout << "bad : " << bad.mu << " " << bad.ed << " "
1069  << bad.pr << " " << bad.en << std::endl;
1070  std::cout << "ret_local,ret: " << ret_local << " "
1071  << ret << std::endl;
1072  std::cout << std::endl;
1073  if (verbose>2) {
1074  char ch;
1075  std::cin >> ch;
1076  }
1077  }
1078 
1079  if (ret_local>ret) {
1080  ret=ret_local;
1081  }
1082 
1083  // End of loop over points in data file
1084  }
1085  // End of temperature loop
1086  }
1087 
1088  dev.mu/=cnt;
1089  dev.ed/=cnt;
1090  dev.pr/=cnt;
1091  dev.en/=cnt;
1092 
1093  if (verbose>0) {
1094  if (k==0) {
1095  std::cout << "Function pair_density(), include rest mass"
1096  << std::endl;
1097  } else if (k==1) {
1098  std::cout << "Function pair_density(), without rest mass"
1099  << std::endl;
1100  } else if (k==2) {
1101  std::cout << "Function pair_density(), include rest mass, "
1102  << "interacting" << std::endl;
1103  } else {
1104  std::cout << "Function pair_density(), without rest mass, "
1105  << "interacting" << std::endl;
1106  }
1107 
1108  std::cout << "Average performance: " << std::endl;
1109  std::cout << "mu: " << dev.mu << " ed: " << dev.ed << " pr: "
1110  << dev.pr << " en: " << dev.en << std::endl;
1111  std::cout << "Worst case: " << std::endl;
1112  std::cout << "mu: " << bad.mu << " ed: " << bad.ed << " pr: "
1113  << bad.pr << " en: " << bad.en << std::endl;
1114  std::cout << "mu: " << mu_bad << " m: " << m_bad
1115  << " T: " << T_bad << " mot: " << mot_bad
1116  << "\n\tpsi: " << psi_bad << std::endl;
1117  std::cout << "ret_local,ret: " << ret_local << " "
1118  << ret << std::endl;
1119  std::cout << std::endl;
1120  if (verbose>2) {
1121  char ch;
1122  std::cin >> ch;
1123  }
1124  }
1125 
1126  if (ret_local>ret) {
1127  ret=ret_local;
1128  }
1129 
1130  // End of k loop
1131  }
1132 
1133  // End of 'if (test_pair)'
1134  }
1135 
1136  // ----------------------------------------------------------------
1137  // Return to the original values
1138 
1139  p=orig;
1140 
1141  return ret;
1142  }
1143 
1144  };
1145 
1146 #ifndef DOXYGEN_NO_O2NS
1147 }
1148 #endif
1149 
1150 #endif
o2scl::part_tl
Particle base class.
Definition: part.h:103
o2scl::table::function_column
void function_column(std::string function, std::string scol)
o2scl::part_calibrate_class::part_calibrate
double part_calibrate(part_t &p, thermo_t &th, bool test_pair, std::string file, bool nr_mode=false, int verbose=0, bool external=false)
Calibrate a particle thermodynamics class with results stored in a table.
Definition: part.h:500
o2scl::table
o2scl::table::get_nlines
size_t get_nlines() const
o2scl::thermo_tl::ed
fp_t ed
Energy density.
Definition: part.h:59
o2scl::part_tl::mu
fp_t mu
Chemical potential.
Definition: part.h:118
o2scl::part_tl::ms
fp_t ms
Effective mass (Dirac unless otherwise specified)
Definition: part.h:122
o2scl::part_tl::type
virtual const char * type()
Return string denoting type ("part_tl")
Definition: part.h:222
o2scl::part_tl::n
fp_t n
Number density.
Definition: part.h:112
o2scl::thermo_tl::thermo_tl
thermo_tl(const thermo_tl &t)
Copy constructor.
Definition: part.h:71
o2scl::part_tl::non_interacting
bool non_interacting
True if the particle is non-interacting (default true)
Definition: part.h:130
o2scl::part_calibrate_class::check_density
void check_density(part1_t &p, part2_t &exact, part3_t &bad, size_t k, double T, double mot, double psi, double &mu_bad, double &m_bad, double &T_bad, double &mot_bad, double &psi_bad, double &ret_local)
Check the density against the exact result and update.
Definition: part.h:304
o2scl::thermo_tl::en
fp_t en
Entropy density.
Definition: part.h:61
o2scl::thermo_tl
A class containing three thermodynamical variables (energy density, pressure, entropy density)
Definition: part.h:52
o2scl::part_tl::en
fp_t en
Entropy density.
Definition: part.h:120
o2scl::part_calibrate_class::set_mass_flags
void set_mass_flags(part_t &p, double mot, double T, size_t k)
Set mass and flags from mot, T, and the index k.
Definition: part.h:248
o2scl_hdf::hdf_file::close
void close()
o2scl::part_tl::operator=
part_tl & operator=(const part_tl &p)
Copy construction with operator=()
Definition: part.h:148
hdf_input
void hdf_input(hdf_file &hf, o2scl::table< vec_t > &t, std::string name)
o2scl::part_tl::g
fp_t g
Degeneracy (e.g. spin and color if applicable)
Definition: part.h:108
o2scl::part_tl::part_tl
part_tl(fp_t mass=0.0, fp_t dof=0.0)
Make a particle of mass mass and degeneracy dof.
Definition: part.h:166
o2scl::part_calibrate_class::check_chem_pot
void check_chem_pot(part1_t &p, part2_t &exact, part3_t &bad, size_t k, double T, double mot, double psi, double &mu_bad, double &m_bad, double &T_bad, double &mot_bad, double &psi_bad, double &ret_local)
Check the chemical potential against the exact result and update.
Definition: part.h:332
o2scl::part_tl::ed
fp_t ed
Energy density.
Definition: part.h:114
o2scl::part_tl::init
virtual void init(fp_t mass, fp_t dof)
Set the mass mass and degeneracy dof.
Definition: part.h:179
o2scl_hdf::hdf_file::open
int open(std::string fname, bool write_access=false, bool err_on_fail=true)
o2scl_settings
lib_settings_class o2scl_settings
o2scl::operator-
thermo operator-(const thermo &left, const thermo &right)
Subtraction operator.
o2scl::part_tl::part_tl
part_tl(const part_tl &p)
Copy constructor.
Definition: part.h:133
o2scl::part_tl::inc_rest_mass
bool inc_rest_mass
If true, include the mass in the energy density and chemical potential (default true)
Definition: part.h:128
o2scl::part_calibrate_class::check_eps
void check_eps(part1_t &p, part2_t &exact, part3_t &bad, size_t k, double T, double mot, double psi, double &mu_bad, double &m_bad, double &T_bad, double &mot_bad, double &psi_bad, double &ret_local)
Check the energy density, pressure, and entropy against the exact result and update.
Definition: part.h:369
o2scl::thermo_tl::pr
fp_t pr
Pressure.
Definition: part.h:57
o2scl::part_tl::m
fp_t m
Mass.
Definition: part.h:110
o2scl::part_tl::pr
fp_t pr
Pressure.
Definition: part.h:116
o2scl::operator+
thermo operator+(const thermo &left, const thermo &right)
Addition operator.
o2scl::thermo_tl::type
const char * type()
Return string denoting type ("thermo")
Definition: part.h:64
o2scl::part_tl::anti
virtual void anti(part_tl &ax)
Make ap an anti-particle with the same mass and degeneracy.
Definition: part.h:205
exc_efilenotfound
exc_efilenotfound
o2scl_hdf::hdf_file
O2SCL_ERR
#define O2SCL_ERR(d, n)
o2scl::thermo_tl::operator=
thermo_tl & operator=(const thermo_tl &t)
Copy construction with operator=()
Definition: part.h:78
o2scl::table::is_column
bool is_column(std::string scol) const
o2scl::part_calibrate_class::set_chem_pot
void set_chem_pot(part_t &p, double psi, double T, size_t k, bool nr_mode)
Set chemical potential from psi, T, the index k, and the flag nr_mode.
Definition: part.h:270
o2scl::part_calibrate_class::check_derivs
void check_derivs(part1_t &p, part2_t &exact, part3_t &bad, size_t k, double T, double mot, double psi, double &mu_bad, double &m_bad, double &T_bad, double &mot_bad, double &psi_bad, double &ret_local)
Check the three derivatives against the exact result and update.
Definition: part.h:429
o2scl::table::get
double get(std::string scol, size_t row) const
o2scl::part_tl::nu
fp_t nu
Effective chemical potential.
Definition: part.h:124
o2scl::part_calibrate_class
Object to organize calibration of particle classes to results stored in a table.
Definition: part.h:241
psi
const double psi

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