Commit 1adff619 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Defined most common physical constants and converted them to the user's favourite system of units.

parent 141b7a94
# Define the system of units to use internally.
UnitSystem:
UnitMass_in_cgs: 1.989e33 # Grams
UnitMass_in_cgs: 1.9885e33 # Grams
UnitLength_in_cgs: 3.0856776e21 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
......
......@@ -28,10 +28,8 @@ import random
# physical constants in cgs
NEWTON_GRAVITY_CGS = 6.672e-8
SOLAR_MASS_IN_CGS = 1.989e33
SOLAR_MASS_IN_CGS = 1.9885e33
PARSEC_IN_CGS = 3.0856776e18
PROTON_MASS_IN_CGS = 1.6726231e24
YEAR_IN_CGS = 3.154e+7
# choice of units
const_unit_length_in_cgs = (1000*PARSEC_IN_CGS)
......
......@@ -280,28 +280,21 @@ int main(int argc, char *argv[]) {
struct UnitSystem us;
struct phys_const prog_const;
units_init(&us, params);
initPhysicalConstants(&us, &prog_const);
init_physical_constants(&us, &prog_const);
if (myrank == 0) {
message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
message("Unit system: U_L = %e cm.", us.UnitLength_in_cgs);
message("Unit system: U_t = %e s.", us.UnitTime_in_cgs);
message("Unit system: U_I = %e A.", us.UnitCurrent_in_cgs);
message("Unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
message("Density units: %e a^%f h^%f.",
units_conversion_factor(&us, UNIT_CONV_DENSITY),
units_a_factor(&us, UNIT_CONV_DENSITY),
units_h_factor(&us, UNIT_CONV_DENSITY));
message("Entropy units: %e a^%f h^%f.",
units_conversion_factor(&us, UNIT_CONV_ENTROPY),
units_a_factor(&us, UNIT_CONV_ENTROPY),
units_h_factor(&us, UNIT_CONV_ENTROPY));
message("Gravity constant = %e", prog_const.newton_gravity);
print_physical_constants(&prog_const);
}
/* Initialise the external potential properties */
struct external_potential potential;
if (with_external_gravity) initPotentialProperties(params, &us, &potential);
if (with_external_gravity && myrank == 0) printPotentialProperties(&potential);
if (with_external_gravity && myrank == 0)
printPotentialProperties(&potential);
/* Read particles and space information from (GADGET) ICs */
char ICfileName[200] = "";
......@@ -390,7 +383,7 @@ int main(int argc, char *argv[]) {
space_map_cells_pre(&s, 0, &map_maxdepth, data);
message("nr of cells at depth %i is %i.", data[0], data[1]);
}
/* Construct the engine policy */
int engine_policies = ENGINE_POLICY | engine_policy_steal;
if (with_hydro) engine_policies |= engine_policy_hydro;
......
......@@ -2602,8 +2602,8 @@ void engine_init(struct engine *e, struct space *s,
e->runners[k].cpuid = k;
e->runners[k].qid = k * nr_queues / e->nr_threads;
}
/* message( "runner %i on cpuid=%i with qid=%i." , e->runners[k].id , */
/* e->runners[k].cpuid , e->runners[k].qid ); */
/* message( "runner %i on cpuid=%i with qid=%i." , e->runners[k].id , */
/* e->runners[k].cpuid , e->runners[k].qid ); */
}
/* Wait for the runner threads to be in place. */
......
......@@ -25,6 +25,7 @@
#include "physical_constants.h"
/* Local headers. */
#include "error.h"
#include "physical_constants_cgs.h"
/**
......@@ -33,10 +34,86 @@
* @param us The current internal system of units.
* @param internal_const The physical constants to initialize.
*/
void initPhysicalConstants(struct UnitSystem* us,
struct phys_const* internal_const) {
void init_physical_constants(struct UnitSystem* us,
struct phys_const* internal_const) {
float dimension[5] = {1, -3, 2, 0, 0};
internal_const->newton_gravity =
NEWTON_GRAVITY_CGS * units_general_conversion_factor(us, dimension);
/* Units are declared as {U_M, U_L, U_t, U_I, U_T} */
const float dimension_G[5] = {-1, 3, -2, 0, 0};
internal_const->const_newton_G =
const_newton_G_cgs / units_general_conversion_factor(us, dimension_G);
const float dimension_c[5] = {0, 1, -1, 0, 0};
internal_const->const_speed_light_c =
const_speed_light_c_cgs /
units_general_conversion_factor(us, dimension_c);
const float dimension_h[5] = {1, -2, -1, 0, 0};
internal_const->const_planck_h =
const_planck_h_cgs / units_general_conversion_factor(us, dimension_h);
internal_const->const_planck_hbar =
const_planck_hbar_cgs / units_general_conversion_factor(us, dimension_h);
const float dimension_k[5] = {1, 2, -2, 0, -1};
internal_const->const_boltzmann_k =
const_boltzmann_k_cgs / units_general_conversion_factor(us, dimension_k);
const float dimension_thomson[5] = {0, 2, 0, 0, 0};
internal_const->const_thomson_cross_section =
const_thomson_cross_section_cgs /
units_general_conversion_factor(us, dimension_thomson);
const float dimension_ev[5] = {1, 2, -2, 0, 0};
internal_const->const_electron_volt =
const_electron_volt_cgs /
units_general_conversion_factor(us, dimension_ev);
const float dimension_charge[5] = {0, 0, -1, 1, 0};
internal_const->const_electron_charge =
const_electron_charge_cgs /
units_general_conversion_factor(us, dimension_charge);
const float dimension_mass[5] = {1, 0, 0, 0, 0};
internal_const->const_electron_mass =
const_electron_mass_cgs /
units_general_conversion_factor(us, dimension_mass);
internal_const->const_proton_mass =
const_proton_mass_cgs /
units_general_conversion_factor(us, dimension_mass);
internal_const->const_solar_mass =
const_solar_mass_cgs /
units_general_conversion_factor(us, dimension_mass);
internal_const->const_earth_mass =
const_earth_mass_cgs /
units_general_conversion_factor(us, dimension_mass);
const float dimension_time[5] = {0, 0, 1, 0, 0};
internal_const->const_year =
const_year_cgs / units_general_conversion_factor(us, dimension_time);
const float dimension_length[5] = {0, 1, 0, 0, 0};
internal_const->const_astronomical_unit =
const_astronomical_unit_cgs /
units_general_conversion_factor(us, dimension_length);
internal_const->const_parsec =
const_parsec_cgs / units_general_conversion_factor(us, dimension_length);
internal_const->const_light_year =
const_light_year_cgs /
units_general_conversion_factor(us, dimension_length);
}
void print_physical_constants(struct phys_const* internal_const) {
message("%25s: %e", "Gravitational constant", internal_const->const_newton_G);
message("%25s: %e", "Speed of light", internal_const->const_speed_light_c);
message("%25s: %e", "Planck constant", internal_const->const_planck_h);
message("%25s: %e", "Boltzmann constant", internal_const->const_boltzmann_k);
message("%25s: %e", "Thomson cross-section",
internal_const->const_thomson_cross_section);
message("%25s: %e", "Electron-Volt", internal_const->const_electron_volt);
message("%25s: %e", "Year", internal_const->const_year);
message("%25s: %e", "Astronomical Unit",
internal_const->const_astronomical_unit);
message("%25s: %e", "Parsec", internal_const->const_parsec);
message("%25s: %e", "Solar mass", internal_const->const_solar_mass);
}
......@@ -28,16 +28,59 @@
/* physical constants in in defined programme units */
struct phys_const {
double newton_gravity;
/* Newton's gravitationl constant */
double const_newton_G;
/* Speed of light in vacuum */
double const_speed_light_c;
/* Planck's constant */
double const_planck_h;
/* Planck's reduced constant */
double const_planck_hbar;
/* Boltzmann's constant */
double const_boltzmann_k;
/* Thomson cross-section */
double const_thomson_cross_section;
/* Charge of the electron */
double const_electron_charge;
/* Electron-Volt */
double const_electron_volt;
/* Mass of the electron */
double const_electron_mass;
/* Mass of the proton */
double const_proton_mass;
/* (Tropical) Year */
double const_year;
/* Astronomical unit */
double const_astronomical_unit;
/* Parsec */
double const_parsec;
/* Light-year */
double const_light_year;
/* Mass of the Sun */
double const_solar_mass;
/* Mass of the Earth */
double const_earth_mass;
};
/**
* @brief Converts physical constants to the internal unit system
*
* @param us The current internal system of units.
* @param internal_const The physical constants to initialize.
*/
void initPhysicalConstants(struct UnitSystem* us,
struct phys_const* internal_const);
void init_physical_constants(struct UnitSystem* us,
struct phys_const* internal_const);
void print_physical_constants(struct phys_const* internal_const);
#endif /* SWIFT_PHYSICAL_CONSTANTS_H */
......@@ -20,14 +20,62 @@
#ifndef SWIFT_PHYSICAL_CONSTANTS_CGS_H
#define SWIFT_PHYSICAL_CONSTANTS_CGS_H
/* physical constants in cgs */
#define NEWTON_GRAVITY_CGS 6.672e-8
#define SOLAR_MASS_IN_CGS 1.989e33
#define PARSEC_IN_CGS 3.0856776e18
#define PROTON_MASS_IN_CGS 1.6726231e24
#define YEAR_IN_CGS 3.154e+7
/* Hydrodynamical constants. */
#define const_hydro_gamma (5.0f / 3.0f)
/* The constants declared in this file should _NOT_ be used directly */
/* Users should use the converted values in the phys_const structure */
/* where all the constants are defined in the system of units specified */
/* in the parameter file. */
/* All values are taken from K.A. Olive et al. (Particle Data Group), Chin. */
/* Phys. C, 38, 090001 (2014) and 2015 update. */
/* http://pdg.lbl.gov/2015/reviews/rpp2015-rev-phys-constants.pdf */
/* http://pdg.lbl.gov/2015/reviews/rpp2015-rev-astrophysical-constants.pdf */
/* Newton's gravitation constant */
const double const_newton_G_cgs = 6.67408e-8; /* g^-1 cm^3 s^-2 */
/* Speed of light in vacuum */
const double const_speed_light_c_cgs = 2.99792458e10; /* cm s^-1 */
/* Planck's constant */
const double const_planck_h_cgs = 6.626070040e-27; /* g cm^-2 s^-1 */
/* Planck's reduced constant */
const double const_planck_hbar_cgs = 1.054571800e-27; /* g cm^-2 s^-1 */
/* Boltzmann's constant */
const double const_boltzmann_k_cgs = 1.38064852e-16; /* g cm^2 s^-2 K^-1 */
/* Thomson cross-section */
const double const_thomson_cross_section_cgs = 6.6524587158e-25; /* cm^2 */
/* Elementary charge */
const double const_electron_charge_cgs = 1.6021766208e-19; /* A s^-1 */
/* Electron-Volt */
const double const_electron_volt_cgs = 1.6021766208e-12; /* g cm^2 s^-2 */
/* Mass of the electron */
const double const_electron_mass_cgs = 9.10938356e-28; /* g */
/* Mass of the proton */
const double const_proton_mass_cgs = 1.672621898e-24; /* g */
/* Tropical year */
const double const_year_cgs = 3.15569252e7; /* s */
/* Astronomical unit */
const double const_astronomical_unit_cgs = 1.49597870700e13; /* cm */
/* Parsec */
const double const_parsec_cgs = 3.08567758149e18; /* cm */
/* Light-year */
const double const_light_year_cgs = 9.46053e17; /* cm */
/* Mass of the Sun */
const double const_solar_mass_cgs = 1.9885e33; /* g */
/* Mass of the Earth */
const double const_earth_mass_cgs = 5.9726e27; /* g */
#endif /* SWIFT_PHYSICAL_CONSTANTS_CGS_H */
......@@ -61,7 +61,7 @@ __attribute__((always_inline))
const struct phys_const* const phys_const,
const struct gpart* const g) {
const float G_newton = phys_const->newton_gravity;
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->point_mass.x;
const float dy = g->x[1] - potential->point_mass.y;
const float dz = g->x[2] - potential->point_mass.z;
......@@ -94,7 +94,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_pointmass(
const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* g) {
const float G_newton = phys_const->newton_gravity;
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->point_mass.x;
const float dy = g->x[1] - potential->point_mass.y;
const float dz = g->x[2] - potential->point_mass.z;
......
......@@ -78,9 +78,10 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
/* #define CELL_ID \ */
/* ((int)(c->loc[0] / DX) * NX *NX + (int)(c->loc[1] / DX) * NX + \ */
/* (int)(c->loc[2] / DX)) */
/* #define OUT \ */
/* if (CELL_ID == MY_CELL) { \ */
/* message(" cell= %d gcount=%d time=%f \n", CELL_ID, c->gcount, r->e->time); \ */
/* #define OUT \ */
/* if (CELL_ID == MY_CELL) { \ */
/* message(" cell= %d gcount=%d time=%f \n", CELL_ID, c->gcount,
* r->e->time); \ */
/* } */
/* #define OUT message(" cell %d %d %f \n",CELL_ID, c->count, r->e->time); */
/* #define OUT if(CELL_ID == MY_CELL) message("\n cell %f %f %f %d %d */
......@@ -165,20 +166,24 @@ void runner_dograv_external(struct runner *r, struct cell *c, int timer) {
/* (g->a_grav[1] * g->a_grav[1]) + */
/* (g->a_grav[2] * g->a_grav[2])) * */
/* dr; */
/* // message("grav_external time= %f\t V_c^2= %f GM/r= %f E= %f L[2]= %f */
/* // x= %f y= %f vx= %f vy= %f\n", r->e->time, v2, fg, E, L[2], g->x[0], */
/* // message("grav_external time= %f\t V_c^2= %f GM/r= %f E= %f L[2]=
* %f */
/* // x= %f y= %f vx= %f vy= %f\n", r->e->time, v2, fg, E, L[2],
* g->x[0], */
/* // g->x[1], vx, vy); */
/* message("%f\t %f %f %f %f %f %f %f %f %f %f %f %f %f\n", r->e->time, */
/* g->tx, g->tv, dt, v2, fg, fga, dr, E, L[2], g->x[0], g->x[1], */
/* message("%f\t %f %f %f %f %f %f %f %f %f %f %f %f %f\n", r->e->time,
*/
/* g->tx, g->tv, dt, v2, fg, fga, dr, E, L[2], g->x[0], g->x[1],
*/
/* vx, vy); */
/* /\* message(" G=%e M=%e\n", r->e->physical_constants->newton_gravity, */
/* /\* message(" G=%e M=%e\n", r->e->physical_constants->newton_gravity,
*/
/* * r->e->potential->point_mass.mass); *\/ */
/* /\* exit(-1); *\/ */
/* } */
}
}
if(timer) TIMER_TOC(timer_dograv_external);
if (timer) TIMER_TOC(timer_dograv_external);
}
/**
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment