diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h index a55fd86ff489fedeb9d2bcb703f62b849ffeeb00..ffa1b9945b5311916995754d8ff385fba524225f 100644 --- a/src/cooling/const_lambda/cooling.h +++ b/src/cooling/const_lambda/cooling.h @@ -35,7 +35,6 @@ #include "physical_constants.h" #include "units.h" - /** * @brief Calculates du/dt in code units for a particle. * @@ -50,13 +49,14 @@ __attribute__((always_inline)) INLINE static float cooling_rate( /* Get particle density */ const float rho = p->rho; - + /* Get cooling function properties */ const float X_H = cooling->hydrogen_mass_abundance; /* Calculate du_dt */ - const float du_dt = -cooling->lambda * ( X_H * rho / phys_const->const_proton_mass) * - ( X_H * rho / phys_const->const_proton_mass) / rho; + const float du_dt = -cooling->lambda * + (X_H * rho / phys_const->const_proton_mass) * + (X_H * rho / phys_const->const_proton_mass) / rho; return du_dt; } @@ -77,7 +77,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( /* Get current internal energy (dt=0) */ const float u_old = hydro_get_internal_energy(p, 0.f); - + /* Internal energy floor */ const float u_floor = cooling->min_energy; @@ -178,9 +178,9 @@ static INLINE void cooling_init_backend( /* convert lambda to code units */ cooling->lambda = lambda_cgs * - units_cgs_conversion_factor(us,UNIT_CONV_TIME) / - (units_cgs_conversion_factor(us,UNIT_CONV_ENERGY) * - units_cgs_conversion_factor(us,UNIT_CONV_VOLUME)); + units_cgs_conversion_factor(us, UNIT_CONV_TIME) / + (units_cgs_conversion_factor(us, UNIT_CONV_ENERGY) * + units_cgs_conversion_factor(us, UNIT_CONV_VOLUME)); } /** @@ -195,8 +195,8 @@ static INLINE void cooling_print_backend( "Cooling function is 'Constant lambda' with " "(lambda,min_energy,hydrogen_mass_abundance,mean_molecular_weight) " "= (%g,%g,%g,%g)", - cooling->lambda, cooling->min_energy, - cooling->hydrogen_mass_abundance, cooling->mean_molecular_weight); + cooling->lambda, cooling->min_energy, cooling->hydrogen_mass_abundance, + cooling->mean_molecular_weight); } #endif /* SWIFT_COOLING_CONST_LAMBDA_H */ diff --git a/src/engine.c b/src/engine.c index b841cdd299dfb27d302c23ba786c43e6d3d26389..e77148ad944cc9c30dfd9e811790f1dff0f2aac9 100644 --- a/src/engine.c +++ b/src/engine.c @@ -131,7 +131,6 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) { const int is_with_cooling = (e->policy & engine_policy_cooling); const int is_with_sourceterms = (e->policy & engine_policy_sourceterms); - /* Are we in a super-cell ? */ if (c->super == c) { @@ -3388,8 +3387,9 @@ void engine_init(struct engine *e, struct space *s, fprintf(e->file_stats, "#%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s " "%14s %14s %14s\n", - "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "E_pot_self", "E_pot_ext", - "E_radcool", "Entropy", "p_x", "p_y", "p_z", "ang_x", "ang_y", "ang_z"); + "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "E_pot_self", + "E_pot_ext", "E_radcool", "Entropy", "p_x", "p_y", "p_z", "ang_x", + "ang_y", "ang_z"); fflush(e->file_stats); char timestepsfileName[200] = ""; diff --git a/src/potential.c b/src/potential.c index 534540ea5bc4e7d05bbd91761548411962f2c327..7c426c18cb1ac5787c563fba5e9722cce9b0c73b 100644 --- a/src/potential.c +++ b/src/potential.c @@ -35,8 +35,7 @@ */ void potential_init(const struct swift_params* parameter_file, const struct phys_const* phys_const, - const struct UnitSystem* us, - const struct space* s, + const struct UnitSystem* us, const struct space* s, struct external_potential* potential) { potential_init_backend(parameter_file, phys_const, us, s, potential); diff --git a/src/potential.h b/src/potential.h index 9b16a0e7aa09a3336d8fa8fe03b77a59ba6cee84..f7fdd0072c502641d36f229337f4616b9033cbda 100644 --- a/src/potential.h +++ b/src/potential.h @@ -48,8 +48,7 @@ /* Now, some generic functions, defined in the source file */ void potential_init(const struct swift_params* parameter_file, const struct phys_const* phys_const, - const struct UnitSystem* us, - const struct space* s, + const struct UnitSystem* us, const struct space* s, struct external_potential* potential); void potential_print(const struct external_potential* potential); diff --git a/src/potential/disc_patch/potential.h b/src/potential/disc_patch/potential.h index 725a7a7634cc34291d6c9d9fdf29acf8c5a035af..130f4feb9aed9c41aee8d6e55c27dbab80c4196b 100644 --- a/src/potential/disc_patch/potential.h +++ b/src/potential/disc_patch/potential.h @@ -33,8 +33,8 @@ #include "parser.h" #include "part.h" #include "physical_constants.h" -#include "units.h" #include "space.h" +#include "units.h" /** * @brief External Potential Properties - Disc patch case @@ -157,16 +157,17 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( * Placeholder for now- just returns 0 * * @param potential The #external_potential used in the run. - * @param phys_const Physical constants in internal units. + * @param phys_const Physical constants in internal units. * @param p Pointer to the particle data. */ - __attribute__((always_inline)) INLINE static float external_gravity_get_potential_energy( +__attribute__((always_inline)) INLINE static float +external_gravity_get_potential_energy( const struct external_potential* potential, const struct phys_const* const phys_const, const struct part* p) { return 0.; - } +} /** * @brief Initialises the external potential properties in the internal system @@ -180,8 +181,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( static INLINE void potential_init_backend( const struct swift_params* parameter_file, const struct phys_const* phys_const, const struct UnitSystem* us, - const struct space* s, - struct external_potential* potential) { + const struct space* s, struct external_potential* potential) { potential->surface_density = parser_get_param_double( parameter_file, "DiscPatchPotential:surface_density"); diff --git a/src/potential/isothermal/potential.h b/src/potential/isothermal/potential.h index bb9206172d245b3bded709b4d04328523d417c47..e77f5d079eb02a2daf9446da6843237b54428759 100644 --- a/src/potential/isothermal/potential.h +++ b/src/potential/isothermal/potential.h @@ -31,8 +31,8 @@ #include "parser.h" #include "part.h" #include "physical_constants.h" -#include "units.h" #include "space.h" +#include "units.h" /** * @brief External Potential Properties - Isothermal sphere case @@ -100,32 +100,33 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep( * @param g Pointer to the g-particle data. */ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( - double time ,const struct external_potential* potential, + double time, const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) { const float dx = g->x[0] - potential->x; const float dy = g->x[1] - potential->y; const float dz = g->x[2] - potential->z; - - const float rinv2 = 1./(dx*dx + dy*dy + dz*dz); + + const float rinv2 = 1. / (dx * dx + dy * dy + dz * dz); const double term = -potential->vrot2_over_G * rinv2; g->a_grav[0] = term * dx; g->a_grav[1] = term * dy; - g->a_grav[2] = term * dz; + g->a_grav[2] = term * dz; } - /** - * @brief Computes the gravitational potential energy of a particle in an isothermal potential. + * @brief Computes the gravitational potential energy of a particle in an + * isothermal potential. * * @param potential The #external_potential used in the run. - * @param phys_const Physical constants in internal units. + * @param phys_const Physical constants in internal units. * @param g Pointer to the particle data. */ - __attribute__((always_inline)) INLINE static float external_gravity_get_potential_energy( +__attribute__((always_inline)) INLINE static float +external_gravity_get_potential_energy( const struct external_potential* potential, const struct phys_const* const phys_const, const struct gpart* g) { @@ -133,8 +134,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( const float dy = g->x[1] - potential->y; const float dz = g->x[2] - potential->z; - return potential->vrot * potential->vrot * 0.5 * log(dx*dx + dy*dy * dz*dz); - } + return potential->vrot * potential->vrot * 0.5 * + log(dx * dx + dy * dy * dz * dz); +} /** * @brief Initialises the external potential properties in the internal system @@ -148,14 +150,16 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( static INLINE void potential_init_backend( const struct swift_params* parameter_file, const struct phys_const* phys_const, const struct UnitSystem* us, - const struct space* s, - struct external_potential* potential) { + const struct space* s, struct external_potential* potential) { - potential->x = s->dim[0]/2. + + potential->x = + s->dim[0] / 2. + parser_get_param_double(parameter_file, "IsothermalPotential:position_x"); - potential->y = s->dim[1]/2. + + potential->y = + s->dim[1] / 2. + parser_get_param_double(parameter_file, "IsothermalPotential:position_y"); - potential->z = s->dim[2]/2. + + potential->z = + s->dim[2] / 2. + parser_get_param_double(parameter_file, "IsothermalPotential:position_z"); potential->vrot = parser_get_param_double(parameter_file, "IsothermalPotential:vrot"); diff --git a/src/potential/none/potential.h b/src/potential/none/potential.h index fb3b27636cd9840f6b0b8fe1b7c6a9e20e851e82..6b9ac7268a396213243abc59fa8e08e097446bfb 100644 --- a/src/potential/none/potential.h +++ b/src/potential/none/potential.h @@ -30,8 +30,8 @@ #include "parser.h" #include "part.h" #include "physical_constants.h" -#include "units.h" #include "space.h" +#include "units.h" /** * @brief External Potential Properties @@ -72,16 +72,17 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( * @brief Computes the gravitational potential energy due to nothing. * * @param potential The #external_potential used in the run. - * @param phys_const Physical constants in internal units. + * @param phys_const Physical constants in internal units. * @param g Pointer to the particle data. */ - __attribute__((always_inline)) INLINE static float external_gravity_get_potential_energy( +__attribute__((always_inline)) INLINE static float +external_gravity_get_potential_energy( const struct external_potential* potential, const struct phys_const* const phys_const, const struct part* g) { return 0.; - } +} /** * @brief Initialises the external potential properties in the internal system @@ -97,8 +98,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( static INLINE void potential_init_backend( const struct swift_params* parameter_file, const struct phys_const* phys_const, const struct UnitSystem* us, - const struct space* s, - struct external_potential* potential) {} + const struct space* s, struct external_potential* potential) {} /** * @brief Prints the properties of the external potential to stdout. diff --git a/src/potential/point_mass/potential.h b/src/potential/point_mass/potential.h index 27ed2bcd4be1d2d85800a9e296e327755e90dc4c..d6926e6dc08a90149a15255f409fdc86d311bb57 100644 --- a/src/potential/point_mass/potential.h +++ b/src/potential/point_mass/potential.h @@ -31,8 +31,8 @@ #include "parser.h" #include "part.h" #include "physical_constants.h" -#include "units.h" #include "space.h" +#include "units.h" /** * @brief External Potential Properties - Point mass case @@ -116,25 +116,26 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( g->a_grav[2] += -potential->mass * dz * rinv3; } - /** - * @brief Computes the gravitational potential energy of a particle in a point mass potential. + * @brief Computes the gravitational potential energy of a particle in a point + * mass potential. * * @param potential The #external_potential used in the run. - * @param phys_const Physical constants in internal units. + * @param phys_const Physical constants in internal units. * @param g Pointer to the particle data. */ - __attribute__((always_inline)) INLINE static float external_gravity_get_potential_energy( +__attribute__((always_inline)) INLINE static float +external_gravity_get_potential_energy( const struct external_potential* potential, const struct phys_const* const phys_const, const struct part* g) { const float dx = g->x[0] - potential->x; const float dy = g->x[1] - potential->y; const float dz = g->x[2] - potential->z; - const float rinv = 1./sqrtf(dx * dx + dy * dy + dz * dz); + const float rinv = 1. / sqrtf(dx * dx + dy * dy + dz * dz); return -phys_const->const_newton_G * potential->mass * r_inv; - } +} /** * @brief Initialises the external potential properties in the internal system @@ -148,8 +149,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( static INLINE void potential_init_backend( const struct swift_params* parameter_file, const struct phys_const* phys_const, const struct UnitSystem* us, - const struct space* s, - struct external_potential* potential) { + const struct space* s, struct external_potential* potential) { potential->x = parser_get_param_double(parameter_file, "PointMassPotential:position_x"); diff --git a/src/potential/softened_isothermal/potential.h b/src/potential/softened_isothermal/potential.h index 9533986f3b6bb6e98f4aba40ae09065e967d71cb..8f7d7cf270353a0a92f5cced8b504a5249272c8f 100644 --- a/src/potential/softened_isothermal/potential.h +++ b/src/potential/softened_isothermal/potential.h @@ -32,8 +32,8 @@ #include "parser.h" #include "part.h" #include "physical_constants.h" -#include "units.h" #include "space.h" +#include "units.h" /** * @brief External Potential Properties - Softened Isothermal sphere case @@ -46,10 +46,12 @@ struct external_potential { /*! Rotation velocity */ double vrot; - /*! Square of vrot, the circular velocity which defines the isothermal potential */ + /*! Square of vrot, the circular velocity which defines the isothermal + * potential */ double vrot2_over_G; - /*! Square of the softening length. Acceleration tends to zero within this distance from the origin */ + /*! Square of the softening length. Acceleration tends to zero within this + * distance from the origin */ double epsilon2; /*! Time-step condition pre-factor */ @@ -74,17 +76,18 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep( const float dy = g->x[1] - potential->y; const float dz = g->x[2] - potential->z; - const float r2_plus_epsilon2_inv = 1.f / (dx * dx + dy * dy + dz * dz + potential->epsilon2); + const float r2_plus_epsilon2_inv = + 1.f / (dx * dx + dy * dy + dz * dz + potential->epsilon2); const float drdv = dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]); const double vrot = potential->vrot; - const float dota_x = - vrot * vrot * r2_plus_epsilon2_inv * (g->v_full[0] - 2.f * drdv * dx * r2_plus_epsilon2_inv); - const float dota_y = - vrot * vrot * r2_plus_epsilon2_inv * (g->v_full[1] - 2.f * drdv * dy * r2_plus_epsilon2_inv); - const float dota_z = - vrot * vrot * r2_plus_epsilon2_inv * (g->v_full[2] - 2.f * drdv * dz * r2_plus_epsilon2_inv); + const float dota_x = vrot * vrot * r2_plus_epsilon2_inv * + (g->v_full[0] - 2.f * drdv * dx * r2_plus_epsilon2_inv); + const float dota_y = vrot * vrot * r2_plus_epsilon2_inv * + (g->v_full[1] - 2.f * drdv * dy * r2_plus_epsilon2_inv); + const float dota_z = vrot * vrot * r2_plus_epsilon2_inv * + (g->v_full[2] - 2.f * drdv * dz * r2_plus_epsilon2_inv); const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z; const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] + g->a_grav[2] * g->a_grav[2]; @@ -97,7 +100,7 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep( * * Note that the accelerations are multiplied by Newton's G constant * later on. - * + * * a = v_rot^2 * (x,y,z) / (r^2 + epsilon^2) * @param time The current time. * @param potential The #external_potential used in the run. @@ -111,24 +114,27 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( const float dx = g->x[0] - potential->x; const float dy = g->x[1] - potential->y; const float dz = g->x[2] - potential->z; - const float r2_plus_epsilon2_inv = 1.f / (dx * dx + dy * dy + dz * dz + potential->epsilon2); + const float r2_plus_epsilon2_inv = + 1.f / (dx * dx + dy * dy + dz * dz + potential->epsilon2); const double term = -potential->vrot2_over_G * r2_plus_epsilon2_inv; g->a_grav[0] = term * dx; g->a_grav[1] = term * dy; - g->a_grav[2] = term * dz; + g->a_grav[2] = term * dz; } /** - * @brief Computes the gravitational potential energy of a particle in an isothermal potential. + * @brief Computes the gravitational potential energy of a particle in an + * isothermal potential. * * @param potential The #external_potential used in the run. - * @param phys_const Physical constants in internal units. + * @param phys_const Physical constants in internal units. * @param g Pointer to the particle data. */ - __attribute__((always_inline)) INLINE static float external_gravity_get_potential_energy( +__attribute__((always_inline)) INLINE static float +external_gravity_get_potential_energy( const struct external_potential* potential, const struct phys_const* const phys_const, const struct part* g) { @@ -136,8 +142,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( const float dy = g->x[1] - potential->y; const float dz = g->x[2] - potential->z; - return potential->vrot * potential->vrot * 0.5 * log(dx*dx + dy*dy * dz*dz + potential->epsilon2); - } + return potential->vrot * potential->vrot * 0.5 * + log(dx * dx + dy * dy * dz * dz + potential->epsilon2); +} /** * @brief Initialises the external potential properties in the internal system * of units. @@ -150,22 +157,25 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( static INLINE void potential_init_backend( const struct swift_params* parameter_file, const struct phys_const* phys_const, const struct UnitSystem* us, - const struct space* s, - struct external_potential* potential) { - - potential->x = s->dim[0]/2. + - parser_get_param_double(parameter_file, "SoftenedIsothermalPotential:position_x"); - potential->y = s->dim[1]/2. + - parser_get_param_double(parameter_file, "SoftenedIsothermalPotential:position_y"); - potential->z = s->dim[2]/2. + - parser_get_param_double(parameter_file, "SoftenedIsothermalPotential:position_z"); - potential->vrot = - parser_get_param_double(parameter_file, "SoftenedIsothermalPotential:vrot"); + const struct space* s, struct external_potential* potential) { + + potential->x = s->dim[0] / 2. + + parser_get_param_double( + parameter_file, "SoftenedIsothermalPotential:position_x"); + potential->y = s->dim[1] / 2. + + parser_get_param_double( + parameter_file, "SoftenedIsothermalPotential:position_y"); + potential->z = s->dim[2] / 2. + + parser_get_param_double( + parameter_file, "SoftenedIsothermalPotential:position_z"); + potential->vrot = parser_get_param_double(parameter_file, + "SoftenedIsothermalPotential:vrot"); potential->timestep_mult = parser_get_param_float( parameter_file, "SoftenedIsothermalPotential:timestep_mult"); const double epsilon = parser_get_param_float( parameter_file, "SoftenedIsothermalPotential:epsilon"); - potential->vrot2_over_G = potential->vrot * potential->vrot / phys_const->const_newton_G; + potential->vrot2_over_G = + potential->vrot * potential->vrot / phys_const->const_newton_G; potential->epsilon2 = epsilon * epsilon; } @@ -182,7 +192,7 @@ static INLINE void potential_print_backend( "%e, %e), vrot = %e " "timestep multiplier = %e, epsilon = %e", potential->x, potential->y, potential->z, potential->vrot, - potential->timestep_mult,sqrtf(potential->epsilon2)); + potential->timestep_mult, sqrtf(potential->epsilon2)); } #endif /* SWIFT_POTENTIAL_ISOTHERMAL_H */ diff --git a/src/space.c b/src/space.c index d2dc6cf511a5c5d8cb9a83dc563f88826836f645..a907f880cd13cb7a4dfacacd029589f61e406009 100644 --- a/src/space.c +++ b/src/space.c @@ -1814,10 +1814,11 @@ void space_init(struct space *s, const struct swift_params *params, } else { for (size_t k = 0; k < Npart; k++) for (int j = 0; j < 3; j++) - if (parts[k].x[j] < 0 || parts[k].x[j] >= dim[j]){ - printf("parts[%lld].x[%d] = %f , dim[%d] = %f\n" , k , j , parts[k].x[j] , j , dim[j]); + if (parts[k].x[j] < 0 || parts[k].x[j] >= dim[j]) { + printf("parts[%lld].x[%d] = %f , dim[%d] = %f\n", k, j, + parts[k].x[j], j, dim[j]); error("Not all particles are within the specified domain."); - } + } } /* Same for the gparts */ diff --git a/src/statistics.c b/src/statistics.c index 96d73b0e610d616dc05aff4dec76961221db1c1f..994b7bb4456d63ed85083c5eb77adc0d35ff065e 100644 --- a/src/statistics.c +++ b/src/statistics.c @@ -123,8 +123,7 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) { const struct part *p = &parts[k]; const struct xpart *xp = &xparts[k]; struct gpart *gp = NULL; - if (p->gpart != NULL) - gp = p->gpart; + if (p->gpart != NULL) gp = p->gpart; /* Get useful variables */ const float dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase; @@ -157,7 +156,8 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) { /* Collect energies. */ stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); stats.E_pot_self += 0.; - stats.E_pot_ext += m * external_gravity_get_potential_energy(potential,phys_const,gp); + stats.E_pot_ext += + m * external_gravity_get_potential_energy(potential, phys_const, gp); stats.E_int += m * hydro_get_internal_energy(p, dt); stats.E_rad += cooling_get_radiated_energy(xp); @@ -230,7 +230,8 @@ void stats_collect_gpart_mapper(void *map_data, int nr_gparts, /* Collect energies. */ stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); stats.E_pot_self += 0.; - stats.E_pot_ext += m * external_gravity_get_potential_energy(potential,phys_const,gp); + stats.E_pot_ext += + m * external_gravity_get_potential_energy(potential, phys_const, gp); } /* Now write back to memory */ @@ -279,9 +280,9 @@ void stats_print_to_file(FILE *file, const struct statistics *stats, " %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e " "%14e %14e %14e\n", time, stats->mass, E_tot, stats->E_kin, stats->E_int, E_pot, - stats->E_pot_self, stats->E_pot_ext, stats->E_rad, stats->entropy, - stats->mom[0], stats->mom[1], stats->mom[2], stats->ang_mom[0], - stats->ang_mom[1], stats->ang_mom[2]); + stats->E_pot_self, stats->E_pot_ext, stats->E_rad, stats->entropy, + stats->mom[0], stats->mom[1], stats->mom[2], stats->ang_mom[0], + stats->ang_mom[1], stats->ang_mom[2]); fflush(file); } diff --git a/src/units.h b/src/units.h index d9b88bfdf22818966b2419d2cae54b78705eb1cf..8761f12d6976125537545e14ad12e88db0f36a6f 100644 --- a/src/units.h +++ b/src/units.h @@ -60,7 +60,6 @@ enum BaseUnits { UNIT_TIME = 2, UNIT_CURRENT = 3, - UNIT_TEMPERATURE = 4 };