diff --git a/examples/UniformBox/makeIC.py b/examples/UniformBox/makeIC.py index 66908cfb02063d175d5980b653c9e706547a7c52..5a907333043951fea0e20748ef6ce4ffc99ec1d7 100644 --- a/examples/UniformBox/makeIC.py +++ b/examples/UniformBox/makeIC.py @@ -29,8 +29,8 @@ from numpy import * periodic= 1 # 1 For periodic box boxSize = 1.0 L = int(sys.argv[1]) # Number of particles along one axis -rho = 1.0e-24 # Density (roughly 1 atom per cubic centimetre) -P = 1.0e-12 # Pressure (at approx 10000K) +rho = 3.3e6 # Density in solar masses per cubic kiloparsec (0.1 hydrogen atoms per cm^3) +P = 1.4e8 # Pressure in code units (at 10^5K) gamma = 5./3. # Gas adiabatic index eta = 1.2349 # 48 ngbs with cubic spline kernel fileName = "uniformBox.hdf5" @@ -38,6 +38,7 @@ fileName = "uniformBox.hdf5" #--------------------------------------------------- numPart = L**3 mass = boxSize**3 * rho / numPart +print mass internalEnergy = P / ((gamma - 1.)*rho) #-------------------------------------------------- diff --git a/examples/UniformBox/uniformBox.yml b/examples/UniformBox/uniformBox.yml index 9f168d39e622327da737c5e76ea88f61e1c13cb6..670bbd483b9fa68f56e26a675df7134114f807b6 100644 --- a/examples/UniformBox/uniformBox.yml +++ b/examples/UniformBox/uniformBox.yml @@ -1,8 +1,8 @@ # Define the system of units to use internally. InternalUnitSystem: - UnitMass_in_cgs: 1 # Grams - UnitLength_in_cgs: 1 # Centimeters - UnitVelocity_in_cgs: 1 # Centimeters per second + UnitMass_in_cgs: 2.0e33 # Solar masses + UnitLength_in_cgs: 3.01e21 # Kilparsecs + UnitVelocity_in_cgs: 1.0e5 # Kilometres per second UnitCurrent_in_cgs: 1 # Amperes UnitTemp_in_cgs: 1 # Kelvin @@ -45,6 +45,8 @@ PointMass: # Cooling parameters (Creasey cooling) (always in cgs) CreaseyCooling: - lambda: 1.0e-22 - min_energy: 0.0 - cooling_tstep_mult: 1.0 \ No newline at end of file + Lambda: 1.0e-22 + minimum_temperature: 1.0e4 + mean_molecular_weight: 0.59 + hydrogen_mass_abundance: 0.75 + cooling_tstep_mult: 1.0 \ No newline at end of file diff --git a/examples/main.c b/examples/main.c index d0f968ec4846322b74638e2cab9e7fae30bd7e8b..85d2f885c0e99c23fee83e2a212f135452c080f2 100644 --- a/examples/main.c +++ b/examples/main.c @@ -342,7 +342,7 @@ int main(int argc, char *argv[]) { /* Initialise the external potential properties */ struct cooling_data cooling; - if (with_cooling) cooling_init(params, &us, &cooling); + if (with_cooling) cooling_init(params, &us, &prog_const, &cooling); if (with_cooling && myrank == 0) cooling_print(&cooling); /* Read particles and space information from (GADGET) ICs */ diff --git a/src/cooling.c b/src/cooling.c index 923af2f44746660ba793c8f09cd5da22ba62dfab..3edb7fcf7af8f4ec83ad16594dc887f7dd176e61 100644 --- a/src/cooling.c +++ b/src/cooling.c @@ -38,6 +38,7 @@ */ void cooling_init(const struct swift_params* parameter_file, struct UnitSystem* us, + const struct phys_const* const phys_const, struct cooling_data* cooling) { #ifdef CONST_COOLING @@ -47,12 +48,21 @@ void cooling_init(const struct swift_params* parameter_file, #endif /* CONST_COOLING */ #ifdef CREASEY_COOLING - cooling->creasey_cooling.lambda = parser_get_param_double(parameter_file, "CreaseyCooling:lambda"); - cooling->creasey_cooling.min_energy = parser_get_param_double(parameter_file, "CreaseyCooling:min_energy"); + cooling->creasey_cooling.lambda = parser_get_param_double(parameter_file, "CreaseyCooling:Lambda"); + cooling->creasey_cooling.min_temperature = parser_get_param_double(parameter_file, "CreaseyCooling:minimum_temperature"); + cooling->creasey_cooling.mean_molecular_weight = parser_get_param_double(parameter_file, "CreaseyCooling:mean_molecular_weight"); + cooling->creasey_cooling.hydrogen_mass_abundance = parser_get_param_double(parameter_file, "CreaseyCooling:hydrogen_mass_abundance"); cooling->creasey_cooling.cooling_tstep_mult = parser_get_param_double(parameter_file, "CreaseyCooling:cooling_tstep_mult"); -#endif /* CREASEY_COOLING */ + /*convert minimum temperature into minimum internal energy*/ + double k_b = phys_const->const_boltzmann_k; + double T_min = cooling->creasey_cooling.min_temperature; + double mu = cooling->creasey_cooling.mean_molecular_weight; + double m_p = phys_const->const_proton_mass; + + cooling->creasey_cooling.min_internal_energy = k_b * T_min / (hydro_gamma_minus_one * mu * m_p); +#endif /* CREASEY_COOLING */ } /** @@ -62,7 +72,7 @@ void cooling_init(const struct swift_params* parameter_file, */ void cooling_print(const struct cooling_data* cooling) { -#ifdef CONST_COOLING */ +#ifdef CONST_COOLING message( "Cooling properties are (lambda, min_energy, tstep multiplier) %g %g %g ", cooling->const_cooling.lambda, @@ -72,16 +82,17 @@ void cooling_print(const struct cooling_data* cooling) { #ifdef CREASEY_COOLING message( - "Cooling properties for Creasey cooling are (lambda, min_energy, tstep multiplier) %g %g %g ", + "Cooling properties for Creasey cooling are (lambda, min_temperature, hydrogen_mass_abundance, mean_molecular_weight, tstep multiplier) %g %g %g %g %g", cooling->creasey_cooling.lambda, - cooling->creasey_cooling.min_energy, + cooling->creasey_cooling.min_temperature, + cooling->creasey_cooling.hydrogen_mass_abundance, + cooling->creasey_cooling.mean_molecular_weight, cooling->creasey_cooling.cooling_tstep_mult); #endif /* CREASEY_COOLING */ } -void update_entropy(const struct cooling_data* cooling, - const struct phys_const* const phys_const, struct part* p, - double dt){ +void update_entropy(const struct phys_const* const phys_const, const struct UnitSystem* us, + const struct cooling_data* cooling, struct part* p, double dt){ /*updates the entropy of a particle after integrating the cooling equation*/ float u_old; @@ -92,7 +103,7 @@ void update_entropy(const struct cooling_data* cooling, // u_old = old_entropy/(GAMMA_MINUS1) * pow(rho,GAMMA_MINUS1); u_old = hydro_get_internal_energy(p,0); // dt = 0 because using current entropy - u_new = calculate_new_thermal_energy(u_old,rho,dt,phys_const,cooling); + u_new = calculate_new_thermal_energy(u_old,rho,dt,cooling,phys_const,us); new_entropy = u_new*pow_minus_gamma_minus_one(rho) * hydro_gamma_minus_one; p->entropy = new_entropy; } @@ -101,8 +112,9 @@ void update_entropy(const struct cooling_data* cooling, thermal energy, density and the timestep dt. Returns the final internal energy*/ float calculate_new_thermal_energy(float u_old, float rho, float dt, - const struct phys_const* const phys_const, - const struct cooling_data* cooling){ + const struct cooling_data* cooling, + const struct phys_const* const phys_const, + const struct UnitSystem* us){ #ifdef CONST_COOLING /*du/dt = -lambda, independent of density*/ float du_dt = -cooling->const_cooling.lambda; @@ -117,19 +129,33 @@ float calculate_new_thermal_energy(float u_old, float rho, float dt, #endif /*CONST_COOLING*/ #ifdef CREASEY_COOLING - /* rho*du/dt = -lambda*n_H^2, where rho = m_p * n_H*/ - float m_p = phys_const->const_proton_mass; - float n_H = rho / m_p; - float lambda = cooling->creasey_cooling.lambda; - float du_dt = -lambda * n_H * n_H / rho; - float u_floor = cooling->creasey_cooling.min_energy; + /* rho*du/dt = -lambda*n_H^2 */ float u_new; - if (u_old - du_dt*dt > u_floor){ - u_new = u_old + du_dt*dt; + float m_p = phys_const->const_proton_mass; + float X_H = cooling->creasey_cooling.hydrogen_mass_abundance; + float lambda = cooling->creasey_cooling.lambda; //this is always in cgs + float u_floor = cooling->creasey_cooling.min_internal_energy; + + /*convert from internal code units to cgs*/ + float dt_cgs = dt * units_cgs_conversion_factor(us,UNIT_CONV_TIME); + float rho_cgs = rho * units_cgs_conversion_factor(us,UNIT_CONV_DENSITY); + float m_p_cgs = m_p * units_cgs_conversion_factor(us,UNIT_CONV_MASS); + float n_H_cgs = rho_cgs / m_p_cgs; + float u_old_cgs = u_old * units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS); + float u_floor_cgs = u_floor * units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS); + float du_dt = -lambda * n_H_cgs * n_H_cgs / rho; + float u_new_cgs; + + if (u_old_cgs - du_dt*dt_cgs > u_floor_cgs){ + u_new_cgs = u_old_cgs + du_dt*dt_cgs; } else{ - u_new = u_floor; + u_new_cgs = u_floor_cgs; } + + /*convert back to internal code units when returning new internal energy*/ + + u_new = u_new_cgs / units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS); #endif /*CREASEY_COOLING*/ return u_new; } diff --git a/src/cooling.h b/src/cooling.h index cbe86bc15c7719910f86f7c11b7854d5462c5f91..e625bbe550275c7e9e2cfd3eb7cbe41dde01317a 100644 --- a/src/cooling.h +++ b/src/cooling.h @@ -52,7 +52,10 @@ struct cooling_data { #ifdef CREASEY_COOLING struct { double lambda; - double min_energy; + double min_temperature; + double hydrogen_mass_abundance; + double mean_molecular_weight; + double min_internal_energy; double cooling_tstep_mult; } creasey_cooling; #endif @@ -84,13 +87,14 @@ cooling_timestep(const struct cooling_data* cooling, /* Now, some generic functions, defined in the source file */ void cooling_init(const struct swift_params* parameter_file, struct UnitSystem* us, + const struct phys_const* const phys_const, struct cooling_data* cooling); void cooling_print(const struct cooling_data* cooling); +void update_entropy(const struct phys_const* const phys_const, const struct UnitSystem* us, + const struct cooling_data* cooling, struct part* p, double dt); float calculate_new_thermal_energy(float u_old, float rho, float dt, - const struct phys_const* const phys_const, - const struct cooling_data* cooling); -void update_entropy(const struct cooling_data* cooling, - const struct phys_const* const phys_const, struct part* p, - double dt); + const struct cooling_data* cooling, + const struct phys_const* const phys_const, + const struct UnitSystem* us); #endif /* SWIFT_COOLING_H */ diff --git a/src/runner.c b/src/runner.c index c53709bf533241c3e32fcf5d42835e16fd8fc87c..d72c47d0c60927c578d8659fe2b715fedc4ddbe9 100644 --- a/src/runner.c +++ b/src/runner.c @@ -151,6 +151,7 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) { const int ti_current = r->e->ti_current; const struct cooling_data *cooling = r->e->cooling_data; const struct phys_const *constants = r->e->physical_constants; + const struct UnitSystem *us = r->e->internalUnits; const double timeBase = r->e->timeBase; double dt; @@ -176,7 +177,7 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) { /* Kick has already updated ti_end, so need to check ti_begin */ if (p->ti_begin == ti_current) { dt = (p->ti_end - p->ti_begin)*timeBase; - update_entropy(cooling, constants, p, dt); + update_entropy(constants, us, cooling ,p, dt); } }