17 #include "gsl/gsl_errno.h"
18 #include "gsl/gsl_integration.h"
19 #include "gsl/gsl_multiroots.h"
31 double degeneracy,
double mass,
double temperature,
32 double effective_chemical_potential,
double statistics,
34 const auto integral = (*integrator)(0.0, 1.0, [&](
double x) {
35 return ((1.0 - x) * (1.0 - x) / (x * x * x * x)) *
39 const double integrated_number_density =
40 (degeneracy / (2.0 * M_PI * M_PI)) * integral;
42 return integrated_number_density;
51 double degeneracy,
double mass,
double number_density,
double temperature,
52 double effective_chemical_potential,
double statistics,
54 const double integrated_number_density =
58 return number_density - integrated_number_density;
62 const gsl_vector *roots_array,
void *parameters, gsl_vector *
function) {
67 const double mass = (par->
mass);
75 gsl_vector_set(
function, 0,
77 degeneracy, mass, number_density, temperature,
84 unsigned int iter, gsl_multiroot_fsolver *solver) {
86 "\n***\nfind_effective_chemical_potential(): iter = %3u \t"
89 iter, gsl_vector_get(solver->x, 0), gsl_vector_get(solver->f, 0));
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;
102 const size_t problem_dimension = 1;
105 degeneracy, mass, number_density, temperature, statistics, integrator};
107 gsl_multiroot_function EffectiveChemicalPotential = {
110 problem_dimension, ¶meters};
112 double roots_array_initial[problem_dimension] = {mu_initial_guess};
114 gsl_vector *roots_array = gsl_vector_alloc(problem_dimension);
115 gsl_vector_set(roots_array, 0, roots_array_initial[0]);
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,
131 status = gsl_multiroot_fsolver_iterate(Root_finder);
139 std::cout <<
"\nGSL error message:\n"
140 << gsl_strerror(status) <<
"\n\n"
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:"
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,
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");
162 status = gsl_multiroot_test_residual(Root_finder->f, solution_precision);
164 if (status == GSL_SUCCESS) {
167 }
while (status == GSL_CONTINUE && iter < 100000);
169 gsl_multiroot_fsolver_free(Root_finder);
170 gsl_vector_free(roots_array);
181 double degeneracy,
double mass,
double number_density,
double temperature,
182 double statistics,
double solution_precision) {
183 assert(temperature > 0);
184 assert(degeneracy > 0);
191 double initial_guess = 0.0;
192 if (statistics == -1.0) {
197 initial_guess = 0.999999 * mass;
199 initial_guess = 1.05 * mass;
203 const double integration_precision = 0.1 * solution_precision;
213 double chemical_potential[1];
214 chemical_potential[0] = 0.0;
220 degeneracy, mass, number_density, temperature, statistics, initial_guess,
221 solution_precision, &
integrator_, chemical_potential);
233 const double precision_fm = 1e-3;
235 double nDensityCheck =
239 if (std::abs((number_density - nDensityCheck) / conversion_factor) >
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 ",
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),
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,
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");
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."
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,
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.");
301 const double small_p = 0.005;
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."
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,
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.");
327 return chemical_potential[0];
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.
void set_precision(double absolute, double relative)
Set precision for the absolute and relative error of the integration.
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.
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.
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