Skip to content
Snippets Groups Projects
Commit f06e74fa authored by Stefan Arridge's avatar Stefan Arridge
Browse files

Improved cooling functions

parent 4bef84ae
No related branches found
No related tags found
1 merge request!242Add the cooling infrastructure and the const_du and const_lambda
...@@ -61,11 +61,12 @@ void cooling_print(const struct cooling_data* cooling) { ...@@ -61,11 +61,12 @@ void cooling_print(const struct cooling_data* cooling) {
cooling->const_cooling.cooling_tstep_mult); cooling->const_cooling.cooling_tstep_mult);
#endif /* CONST_COOLING */ #endif /* CONST_COOLING */
} }
int update_entropy(const struct cooling_data* cooling,
const struct phys_const* const phys_const, struct part* p, float dt){ void update_entropy(const struct cooling_data* cooling,
const struct phys_const* const phys_const, struct part* p,
float dt){
/*updates the entropy of a particle after integrating the cooling equation*/ /*updates the entropy of a particle after integrating the cooling equation*/
int status == 0;
float u_old; float u_old;
float u_new; float u_new;
float new_entropy; float new_entropy;
...@@ -73,34 +74,32 @@ int update_entropy(const struct cooling_data* cooling, ...@@ -73,34 +74,32 @@ int update_entropy(const struct cooling_data* cooling,
float rho = p->rho; float rho = p->rho;
u_old = old_entropy/(GAMMA_MINUS1) * pow(rho,GAMMA_MINUS1); u_old = old_entropy/(GAMMA_MINUS1) * pow(rho,GAMMA_MINUS1);
status = calculate_new_thermal_energy(u_old,&u_new,dt,cooling): u_new = calculate_new_thermal_energy(u_old,dt,cooling):
new_entropy = u_new/pow(rho,GAMMA_MINUS1) * GAMMA_MINUS1;
if (status == 0){ p->entropy = new_entropy
new_entropy = u_new/pow(rho,GAMMA_MINUS1) * GAMMA_MINUS1;
p->entropy = new_entropy
}
else
message("Error with cooling, particle's entropy has not been updated");
return status;
} }
#ifdef CONST_COOLING
int calculate_new_thermal_energy(float u_old, float* u_new, float dt, const struct cooling_data* cooling){
float calculate_new_thermal_energy(float u_old, float dt, const struct cooling_data* cooling){
#ifdef CONST_COOLING
//This function integrates the cooling equation, given the initial thermal energy and the timestep dt. //This function integrates the cooling equation, given the initial thermal energy and the timestep dt.
//Returns 0 if successful and 1 if not //Returns 0 if successful and 1 if not
int status = 0; int status = 0;
float du_dt = cooling->const_cooling.lambda; float du_dt = cooling->const_cooling.lambda;
float u_floor = cooling->const_cooling.min_energy; float u_floor = cooling->const_cooling.min_energy;
if (u_old - du_dt*dt > min_energy): float u_new;
*u_new = u_old - du_dt*dt; if (u_old - du_dt*dt > min_energy){
else: u_new = u_old - du_dt*dt;
*u_new = min_energy }
else{
u_new = min_energy;
}
return u_new;
#endif /*CONST_COOLING*/
}
return status}
#endif /*CONST_COOLING
...@@ -61,7 +61,8 @@ struct cooling_data { ...@@ -61,7 +61,8 @@ struct cooling_data {
*/ */
__attribute__((always_inline)) INLINE static float __attribute__((always_inline)) INLINE static float
cooling_timestep(const struct cooling_data* cooling, cooling_timestep(const struct cooling_data* cooling,
const struct phys_const* const phys_const, const struct part* const p) { const struct phys_const* const phys_const,
const struct part* const p) {
const float cooling_rate = get_cooling_rate( p->density, p->internal_energy, cooling ); const float cooling_rate = get_cooling_rate( p->density, p->internal_energy, cooling );
return cooling->const_cooling.cooling_tstep_mult * p->internal_energy / cooling_rate; return cooling->const_cooling.cooling_tstep_mult * p->internal_energy / cooling_rate;
...@@ -93,5 +94,8 @@ void cooling_init(const struct swift_params* parameter_file, ...@@ -93,5 +94,8 @@ void cooling_init(const struct swift_params* parameter_file,
struct cooling_data* cooling); struct cooling_data* cooling);
void cooling_print(const struct cooling_data* cooling); void cooling_print(const struct cooling_data* cooling);
float calculate_new_thermal_energy(float u_old, float dt, const struct cooling_data* cooling);
void update_entropy(const struct cooling_data* cooling,
const struct phys_const* const phys_const, struct part* p,
float dt);
#endif /* SWIFT_COOLING_H */ #endif /* SWIFT_COOLING_H */
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment