Version: SMASH-3.1
chemicalpotential.cc
Go to the documentation of this file.
1 /*
2  *
3  * Copyright (c) 2020-2022
4  * SMASH Team
5  *
6  * GNU General Public License (GPLv3 or later)
7  *
8  */
9 
11 
12 #include <cmath>
13 #include <cstdlib>
14 #include <iostream>
15 #include <sstream>
16 
17 #include "gsl/gsl_errno.h"
18 #include "gsl/gsl_integration.h"
19 #include "gsl/gsl_multiroots.h"
20 
21 #include "smash/constants.h"
22 #include "smash/distributions.h"
23 #include "smash/logging.h"
24 
25 namespace smash {
26 
31  double degeneracy, double mass, double temperature,
32  double effective_chemical_potential, double statistics,
33  Integrator *integrator) {
34  const auto integral = (*integrator)(0.0, 1.0, [&](double x) {
35  return ((1.0 - x) * (1.0 - x) / (x * x * x * x)) *
36  juttner_distribution_func((1.0 - x) / x, mass, temperature,
37  effective_chemical_potential, statistics);
38  });
39  const double integrated_number_density =
40  (degeneracy / (2.0 * M_PI * M_PI)) * integral;
41 
42  return integrated_number_density;
43 }
44 
45 /*
46  * Root equations and procedure for finding the effective chemical potential
47  * for a given particle species.
48  */
49 
51  double degeneracy, double mass, double number_density, double temperature,
52  double effective_chemical_potential, double statistics,
53  Integrator *integrator) {
54  const double integrated_number_density =
55  density_one_species(degeneracy, mass, temperature,
56  effective_chemical_potential, statistics, integrator);
57 
58  return number_density - integrated_number_density;
59 }
60 
62  const gsl_vector *roots_array, void *parameters, gsl_vector *function) {
63  struct ParametersForChemPotSolver *par =
64  static_cast<ParametersForChemPotSolver *>(parameters);
65 
66  const double degeneracy = (par->degeneracy);
67  const double mass = (par->mass);
68  const double number_density = (par->number_density);
69  const double temperature = (par->temperature);
70  const double statistics = (par->statistics);
71  Integrator *integrator = (par->integrator);
72 
73  const double effective_chemical_potential = gsl_vector_get(roots_array, 0);
74 
75  gsl_vector_set(function, 0,
77  degeneracy, mass, number_density, temperature,
78  effective_chemical_potential, statistics, integrator));
79 
80  return GSL_SUCCESS;
81 }
82 
84  unsigned int iter, gsl_multiroot_fsolver *solver) {
85  std::printf(
86  "\n***\nfind_effective_chemical_potential(): iter = %3u \t"
87  "x = % .3f \t"
88  "f(x) = % .3e \n",
89  iter, gsl_vector_get(solver->x, 0), gsl_vector_get(solver->f, 0));
90 }
91 
93  double degeneracy, double mass, double number_density, double temperature,
94  double statistics, double mu_initial_guess, double solution_precision,
95  Integrator *integrator, double *effective_chemical_potential) {
96  const gsl_multiroot_fsolver_type *Solver_name;
97  gsl_multiroot_fsolver *Root_finder;
98 
99  int status;
100  size_t iter = 0;
101 
102  const size_t problem_dimension = 1;
103 
104  struct ParametersForChemPotSolver parameters = {
105  degeneracy, mass, number_density, temperature, statistics, integrator};
106 
107  gsl_multiroot_function EffectiveChemicalPotential = {
110  problem_dimension, &parameters};
111 
112  double roots_array_initial[problem_dimension] = {mu_initial_guess};
113 
114  gsl_vector *roots_array = gsl_vector_alloc(problem_dimension);
115  gsl_vector_set(roots_array, 0, roots_array_initial[0]);
116 
117  Solver_name = gsl_multiroot_fsolver_hybrids;
118  Root_finder = gsl_multiroot_fsolver_alloc(Solver_name, problem_dimension);
119  gsl_multiroot_fsolver_set(Root_finder, &EffectiveChemicalPotential,
120  roots_array);
121 
122  // print_state_effective_chemical_potential (iter, Root_finder);
123 
124  do {
125  iter++;
126 
127  /*
128  * gsl_multiroot_fsolver_iterate returns either 0 for a correct behavior,
129  * or an error code (a positive integer) when the solver is stuck
130  */
131  status = gsl_multiroot_fsolver_iterate(Root_finder);
132 
133  // print_state_effective_chemical_potential (iter, Root_finder);
134 
135  /*
136  * Check whether the solver is stuck
137  */
138  if (status) {
139  std::cout << "\nGSL error message:\n"
140  << gsl_strerror(status) << "\n\n"
141  << std::endl;
142  const double conversion_factor =
144  logg[LogArea::Main::id].warn(
145  "\n\nThe GSL solver\nfind_effective_chemical_potential\nis stuck!"
146  "\n\nInput parameters:"
147  "\n degeneracy = ",
148  degeneracy, "\n mass [GeV] = ", mass,
149  "\nnumber_density [fm^-3] = ", number_density / (conversion_factor),
150  "\n temperature [GeV] = ", temperature,
151  "\n statistics = ", statistics,
152  "\n solution_precision = ", solution_precision,
153  "\n\n"
154  "Initialization cannot continue without calculating the "
155  "chemical potential. \nTry adjusting the initial "
156  "guess or the solution precision.\n\n\n");
157  throw std::runtime_error(
158  "Chemical potential calculation returned no result.\n\n");
159  continue;
160  }
161 
162  status = gsl_multiroot_test_residual(Root_finder->f, solution_precision);
163 
164  if (status == GSL_SUCCESS) {
165  effective_chemical_potential[0] = gsl_vector_get(Root_finder->x, 0);
166  }
167  } while (status == GSL_CONTINUE && iter < 100000);
168 
169  gsl_multiroot_fsolver_free(Root_finder);
170  gsl_vector_free(roots_array);
171 
172  return 0;
173 }
174 
175 /*
176  * A convenience wrapper for finding the effective chemical potential for a
177  * given particle species and performing sanity checks.
178  */
179 
181  double degeneracy, double mass, double number_density, double temperature,
182  double statistics, double solution_precision) {
183  assert(temperature > 0); // zero temperature case is not considered here
184  assert(degeneracy > 0);
185  assert(mass > 0);
186  /*
187  * The initial guess for the GSL chemical potential finding
188  * procedure; the guess has to be different for Bose and Fermi
189  * particles.
190  */
191  double initial_guess = 0.0;
192  if (statistics == -1.0) {
193  /*
194  * Such initial guess allows one to sample really close to the Bose
195  * condensate region; at the same time, it works fine in other regions.
196  */
197  initial_guess = 0.999999 * mass;
198  } else {
199  initial_guess = 1.05 * mass;
200  }
201 
203  const double integration_precision = 0.1 * solution_precision;
204  integrator_.set_precision(integration_precision, integration_precision);
205 
206  /*
207  * Array for the GSL chemical potential finder.
208  * Use of array is motivated by anticipating a generalization or overloading
209  * of this solver to include effective mass calculation, which will lead to
210  * multidimensional root finding, in which case the result is best passed as
211  * an array.
212  */
213  double chemical_potential[1];
214  chemical_potential[0] = 0.0;
215 
216  /*
217  * Call the GSL effective chemical potential finding procedure.
218  */
220  degeneracy, mass, number_density, temperature, statistics, initial_guess,
221  solution_precision, &integrator_, chemical_potential);
222 
223  /*
224  * Sanity checks are performed for the obtained value of the chemical
225  * potential .
226  */
227 
228  /*
229  * Check that the calculated chemical potential reproduces the input number
230  * density within acceptable error.
231  */
232  // precision for number density in units of fm^-3
233  const double precision_fm = 1e-3;
234  const double conversion_factor = smash::hbarc * smash::hbarc * smash::hbarc;
235  double nDensityCheck =
236  density_one_species(degeneracy, mass, temperature, chemical_potential[0],
237  statistics, &integrator_);
238 
239  if (std::abs((number_density - nDensityCheck) / conversion_factor) >
240  precision_fm) {
241  logg[LogArea::Main::id].warn(
242  "\n\nThe calculated chemical potential = ", chemical_potential[0],
243  " [GeV] does not reproduce the input number density to within ",
244  precision_fm,
245  "!"
246  "\nnumber_density [GeV^3] = ",
247  number_density, "\n nDensityCheck [GeV^3] = ", nDensityCheck,
248  "\nnumber_density [fm^-3] = ", number_density / (conversion_factor),
249  "\n nDensityCheck [fm^-3] = ", nDensityCheck / (conversion_factor),
250  "\n\n"
251  "Input parameters:"
252  "\n degeneracy = ",
253  degeneracy, "\n mass [GeV] = ", mass,
254  "\nnumber_density [fm^-3] = ", number_density / (conversion_factor),
255  "\n temperature [GeV] = ", temperature,
256  "\n statistics = ", statistics,
257  "\n solution_precision = ", solution_precision,
258  "\n\n"
259  "Initialization cannot continue without a correct "
260  "calculation of the chemical potential."
261  "\nTry adjusting the initial guess or the solution precision.\n\n\n");
262  throw std::runtime_error(
263  "Chemical potential calculation does not reproduce the "
264  "input number density.\n\n");
265  }
266 
267  /*
268  * Check if the obtained value of the effective chemical potential indicates
269  * that the number density is such that for bosons one has encountered the
270  * Bose-Einstein condensate. We need mu < mass in order for the system NOT to
271  * be in the Bose-Einstein condensate region.
272  */
273  if ((statistics == -1) && (chemical_potential[0] > mass)) {
274  logg[LogArea::Main::id].warn(
275  "\n\nThe calculated chemical potential indicates the "
276  "input number density is such that Bose-Einstein "
277  "condensate is encountered."
278  "Input parameters:"
279  "\n degeneracy = ",
280  degeneracy, "\n mass [GeV] = ", mass,
281  "\nnumber_density [fm^-3] = ", number_density / (conversion_factor),
282  "\n temperature [GeV] = ", temperature,
283  "\n statistics = ", statistics,
284  "\n solution_precision = ", solution_precision,
285  "\n\n"
286  "Initialization cannot proceed with a pathological "
287  "distribution function leading to the Bose-Einstein "
288  "condensate. \nTry increasing the temperature or "
289  "decreasing the number density.\n\n\n");
290  throw std::runtime_error(
291  "Chemical potential calculation indicates that the Bose-Einstein "
292  "condensate is produced.");
293  }
294 
295  /*
296  * Check if the distribution is pathological, i.e., has negative values
297  * for low momenta (Bose statistics, in the Bose-Einstein condensate regime,
298  * may produce pathological distributions for which f(p) < 0 ).
299  */
301  const double small_p = 0.005;
302  const double distribution_at_small_p = juttner_distribution_func(
303  small_p, mass, temperature, chemical_potential[0], statistics);
304  if (distribution_at_small_p < 0.0) {
305  logg[LogArea::Main::id].warn(
306  "\n\nNegative values of the distribution function at "
307  "small momenta indicate that Bose-Einstein "
308  "condensate is encountered."
309  "Input parameters:"
310  "\n degeneracy = ",
311  degeneracy, "\n mass [GeV] = ", mass,
312  "\nnumber_density [fm^-3] = ", number_density / (conversion_factor),
313  "\n temperature [GeV] = ", temperature,
314  "\n statistics = ", statistics,
315  "\n solution_precision = ", solution_precision, "\n\nf(p=", small_p,
316  ") = ", distribution_at_small_p,
317  "\n\n"
318  "Initialization cannot proceed with a pathological "
319  "distribution function leading to the Bose-Einstein "
320  "condensate. \nTry increasing the temperature or "
321  "decreasing the number density.\n\n\n");
322  throw std::runtime_error(
323  "Negative values of the distribution function at low momenta "
324  "indicate that the Bose-Einstein condensate is produced.");
325  }
326 
327  return chemical_potential[0];
328 }
329 
330 } // namespace smash
static double root_equation_effective_chemical_potential(double degeneracy, double mass, double number_density, double temperature, double effective_chemical_potential, double statistics, Integrator *integrator)
Root equation for finding the value of the effective chemical potential for one particle species.
Integrator integrator_
A wrapper for gsl numerical integration.
double effective_chemical_potential(double degeneracy, double mass, double number_density, double temperature, double statistics, double solution_precision)
Convenience wrapper for finding the effective chemical potential for a given particle species and per...
static void print_state_effective_chemical_potential(unsigned int iter, gsl_multiroot_fsolver *solver)
A GSL utility which allows for printing out the status of the solver during the root finding procedur...
static double density_one_species(double degeneracy, double mass, double temperature, double effective_chemical_potential, double statistics, Integrator *integrator)
Vector number density of one particle species, obtained through integrating the Juttner distribution ...
static int root_equation_effective_chemical_potential_for_GSL(const gsl_vector *roots_array, void *parameters, gsl_vector *function)
Root equation for finding the value of the effective chemical potential for one particle species,...
static int find_effective_chemical_potential(double degeneracy, double mass, double number_density, double temperature, double statistics, double mu_initial_guess, double solution_precision, Integrator *integrator, double *effective_chemical_potential)
A GSL root solver for finding the effective chemical potential.
A C++ interface for numerical integration in one dimension with the GSL CQUAD integration functions.
Definition: integrate.h:106
void set_precision(double absolute, double relative)
Set precision for the absolute and relative error of the integration.
Definition: integrate.h:168
Collection of useful constants that are known at compile time.
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > logg
An array that stores all pre-configured Logger objects.
Definition: logging.cc:39
Definition: action.h:24
double juttner_distribution_func(double momentum_radial, double mass, double temperature, double effective_chemical_potential, double statistics)
Relativistic Juttner distribution function is just a convenience wrapper for displaying Fermi,...
constexpr double hbarc
GeV <-> fm conversion factor.
Definition: constants.h:25
Struct, root equations, and procedure for finding the effective chemical potential for a given partic...
double temperature
temperature T of the system [GeV]
double degeneracy
degeneracy g of the particle species
double number_density
number density n of the particle species [GeV^3]
double statistics
statistics quantum statistics of the particles species (+1 for Fermi, -1 for Bose,...
Integrator * integrator
wrapper for gsl numerical integration
double mass
(pole) mass m of the particle species