Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
SWIFT
SWIFTsim
Commits
47883b9f
Commit
47883b9f
authored
Aug 29, 2016
by
Matthieu Schaller
Browse files
Code formatting
parent
12e478a5
Changes
6
Hide whitespace changes
Inline
Side-by-side
src/cooling.c
View file @
47883b9f
...
...
@@ -25,9 +25,9 @@
#include
"../config.h"
/* This object's header. */
#include
"adiabatic_index.h"
#include
"cooling.h"
#include
"hydro.h"
#include
"adiabatic_index.h"
/**
* @brief Initialises the cooling properties in the internal system
...
...
@@ -38,29 +38,40 @@
* @param cooling The cooling properties to initialize
*/
void
cooling_init
(
const
struct
swift_params
*
parameter_file
,
struct
UnitSystem
*
us
,
const
struct
phys_const
*
const
phys_const
,
struct
cooling_data
*
cooling
)
{
struct
UnitSystem
*
us
,
const
struct
phys_const
*
const
phys_const
,
struct
cooling_data
*
cooling
)
{
#ifdef CONST_COOLING
cooling
->
const_cooling
.
lambda
=
parser_get_param_double
(
parameter_file
,
"Cooling:lambda"
);
cooling
->
const_cooling
.
min_energy
=
parser_get_param_double
(
parameter_file
,
"Cooling:min_energy"
);
cooling
->
const_cooling
.
cooling_tstep_mult
=
parser_get_param_double
(
parameter_file
,
"Cooling:cooling_tstep_mult"
);
cooling
->
const_cooling
.
lambda
=
parser_get_param_double
(
parameter_file
,
"Cooling:lambda"
);
cooling
->
const_cooling
.
min_energy
=
parser_get_param_double
(
parameter_file
,
"Cooling:min_energy"
);
cooling
->
const_cooling
.
cooling_tstep_mult
=
parser_get_param_double
(
parameter_file
,
"Cooling:cooling_tstep_mult"
);
#endif
/* CONST_COOLING */
#ifdef CREASEY_COOLING
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"
);
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"
);
/*convert minimum temperature into minimum internal energy*/
float
u_floor
=
phys_const
->
const_boltzmann_k
*
cooling
->
creasey_cooling
.
min_temperature
/
(
hydro_gamma_minus_one
*
cooling
->
creasey_cooling
.
mean_molecular_weight
*
phys_const
->
const_proton_mass
);
float
u_floor_cgs
=
u_floor
*
units_cgs_conversion_factor
(
us
,
UNIT_CONV_ENERGY_PER_UNIT_MASS
);
cooling
->
creasey_cooling
.
min_internal_energy
=
u_floor
;
float
u_floor
=
phys_const
->
const_boltzmann_k
*
cooling
->
creasey_cooling
.
min_temperature
/
(
hydro_gamma_minus_one
*
cooling
->
creasey_cooling
.
mean_molecular_weight
*
phys_const
->
const_proton_mass
);
float
u_floor_cgs
=
u_floor
*
units_cgs_conversion_factor
(
us
,
UNIT_CONV_ENERGY_PER_UNIT_MASS
);
cooling
->
creasey_cooling
.
min_internal_energy
=
u_floor
;
cooling
->
creasey_cooling
.
min_internal_energy_cgs
=
u_floor_cgs
;
#endif
/* CREASEY_COOLING */
}
...
...
@@ -72,58 +83,61 @@ 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
,
cooling
->
const_cooling
.
min_energy
,
cooling
->
const_cooling
.
lambda
,
cooling
->
const_cooling
.
min_energy
,
cooling
->
const_cooling
.
cooling_tstep_mult
);
#endif
/* CONST_COOLING */
#ifdef CREASEY_COOLING
message
(
"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_temperature
,
"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_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
phys_const
*
const
phys_const
,
const
struct
UnitSystem
*
us
,
const
struct
cooling_data
*
cooling
,
struct
part
*
p
,
float
dt
){
void
update_entropy
(
const
struct
phys_const
*
const
phys_const
,
const
struct
UnitSystem
*
us
,
const
struct
cooling_data
*
cooling
,
struct
part
*
p
,
float
dt
)
{
/*updates the entropy of a particle after integrating the cooling equation*/
float
u_old
;
float
u_new
;
float
new_entropy
;
//float old_entropy = p->entropy;
//
float old_entropy = p->entropy;
float
rho
=
p
->
rho
;
// 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
,
cooling
,
phys_const
,
us
);
new_entropy
=
u_new
*
pow_minus_gamma_minus_one
(
rho
)
*
hydro_gamma_minus_one
;
u_old
=
hydro_get_internal_energy
(
p
,
0
);
// dt = 0 because using current entropy
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
;
}
/*This function integrates the cooling equation, given the initial
thermal energy, density and the timestep dt. Returns the final internal energy*/
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
cooling_data
*
cooling
,
const
struct
phys_const
*
const
phys_const
,
const
struct
UnitSystem
*
us
){
float
calculate_new_thermal_energy
(
float
u_old
,
float
rho
,
float
dt
,
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
;
float
u_floor
=
cooling
->
const_cooling
.
min_energy
;
float
u_new
;
if
(
u_old
-
du_dt
*
dt
>
u_floor
){
u_new
=
u_old
+
du_dt
*
dt
;
}
else
{
if
(
u_old
-
du_dt
*
dt
>
u_floor
)
{
u_new
=
u_old
+
du_dt
*
dt
;
}
else
{
u_new
=
u_floor
;
}
#endif
/*CONST_COOLING*/
...
...
@@ -132,33 +146,30 @@ float calculate_new_thermal_energy(float u_old, float rho, float dt,
/* rho*du/dt = -lambda*n_H^2 */
float
u_new
;
float
X_H
=
cooling
->
creasey_cooling
.
hydrogen_mass_abundance
;
float
lambda_cgs
=
cooling
->
creasey_cooling
.
lambda
;
//this is always in cgs
float
lambda_cgs
=
cooling
->
creasey_cooling
.
lambda
;
//
this is always in cgs
float
u_floor_cgs
=
cooling
->
creasey_cooling
.
min_internal_energy_cgs
;
/*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
=
phys_const
->
const_proton_mass
*
units_cgs_conversion_factor
(
us
,
UNIT_CONV_MASS
);
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
=
phys_const
->
const_proton_mass
*
units_cgs_conversion_factor
(
us
,
UNIT_CONV_MASS
);
float
n_H_cgs
=
X_H
*
rho_cgs
/
m_p_cgs
;
float
u_old_cgs
=
u_old
*
units_cgs_conversion_factor
(
us
,
UNIT_CONV_ENERGY_PER_UNIT_MASS
);
float
u_old_cgs
=
u_old
*
units_cgs_conversion_factor
(
us
,
UNIT_CONV_ENERGY_PER_UNIT_MASS
);
float
du_dt_cgs
=
-
lambda_cgs
*
n_H_cgs
*
n_H_cgs
/
rho_cgs
;
float
u_new_cgs
;
if
(
u_old_cgs
+
du_dt_cgs
*
dt_cgs
>
u_floor_cgs
){
u_new_cgs
=
u_old_cgs
+
du_dt_cgs
*
dt_cgs
;
}
else
{
if
(
u_old_cgs
+
du_dt_cgs
*
dt_cgs
>
u_floor_cgs
)
{
u_new_cgs
=
u_old_cgs
+
du_dt_cgs
*
dt_cgs
;
}
else
{
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
);
u_new
=
u_new_cgs
/
units_cgs_conversion_factor
(
us
,
UNIT_CONV_ENERGY_PER_UNIT_MASS
);
#endif
/*CREASEY_COOLING*/
return
u_new
;
}
src/cooling.h
View file @
47883b9f
...
...
@@ -72,30 +72,32 @@ struct cooling_data {
* @param phys_const The physical constants in internal units.
* @param Pointer to the particle data.
*/
__attribute__
((
always_inline
))
INLINE
static
double
cooling_timestep
(
const
struct
cooling_data
*
cooling
,
const
struct
phys_const
*
const
phys_const
,
const
struct
part
*
const
p
)
{
__attribute__
((
always_inline
))
INLINE
static
double
cooling_timestep
(
const
struct
cooling_data
*
cooling
,
const
struct
phys_const
*
const
phys_const
,
const
struct
part
*
const
p
)
{
const
double
cooling_rate
=
cooling
->
const_cooling
.
lambda
;
const
double
internal_energy
=
hydro_get_internal_energy
(
p
,
0
);
// dt = 0 because using current entropy
return
cooling
->
const_cooling
.
cooling_tstep_mult
*
internal_energy
/
cooling_rate
;
const
double
internal_energy
=
hydro_get_internal_energy
(
p
,
0
);
// dt = 0 because using current entropy
return
cooling
->
const_cooling
.
cooling_tstep_mult
*
internal_energy
/
cooling_rate
;
}
#endif
/* CONST_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
);
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
,
float
dt
);
float
calculate_new_thermal_energy
(
float
u_old
,
float
rho
,
float
dt
,
const
struct
cooling_data
*
cooling
,
const
struct
phys_const
*
const
phys_const
,
const
struct
UnitSystem
*
us
);
void
update_entropy
(
const
struct
phys_const
*
const
phys_const
,
const
struct
UnitSystem
*
us
,
const
struct
cooling_data
*
cooling
,
struct
part
*
p
,
float
dt
);
float
calculate_new_thermal_energy
(
float
u_old
,
float
rho
,
float
dt
,
const
struct
cooling_data
*
cooling
,
const
struct
phys_const
*
const
phys_const
,
const
struct
UnitSystem
*
us
);
#endif
/* SWIFT_COOLING_H */
src/engine.c
View file @
47883b9f
...
...
@@ -82,7 +82,7 @@ const char *engine_policy_names[15] = {"none",
"external_gravity"
,
"cosmology_integration"
,
"drift_all"
,
"cooling"
};
"cooling"
};
/** The rank of the engine as a global variable (for messages). */
int
engine_rank
;
...
...
@@ -186,8 +186,7 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
struct
scheduler
*
s
=
&
e
->
sched
;
const
int
is_fixdt
=
(
e
->
policy
&
engine_policy_fixdt
)
==
engine_policy_fixdt
;
const
int
is_with_cooling
=
(
e
->
policy
&
engine_policy_cooling
)
==
engine_policy_cooling
;
(
e
->
policy
&
engine_policy_cooling
)
==
engine_policy_cooling
;
/* Is this the super-cell? */
if
(
super
==
NULL
&&
(
c
->
density
!=
NULL
||
(
c
->
count
>
0
&&
!
c
->
split
)))
{
...
...
@@ -225,9 +224,8 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
#endif
if
(
is_with_cooling
)
c
->
cooling
=
scheduler_addtask
(
s
,
task_type_cooling
,
task_subtype_none
,
0
,
0
,
c
,
NULL
,
0
);
c
->
cooling
=
scheduler_addtask
(
s
,
task_type_cooling
,
task_subtype_none
,
0
,
0
,
c
,
NULL
,
0
);
}
}
...
...
@@ -1793,10 +1791,10 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
scheduler_addunlock
(
sched
,
t
->
ci
->
init
,
t
);
scheduler_addunlock
(
sched
,
t
,
t
->
ci
->
kick
);
}
/* Cooling tasks should depend on kick and does not unlock anything since
it is the last task*/
else
if
(
t
->
type
==
task_type_cooling
)
{
else
if
(
t
->
type
==
task_type_cooling
)
{
scheduler_addunlock
(
sched
,
t
->
ci
->
kick
,
t
);
}
}
...
...
@@ -2780,7 +2778,7 @@ void engine_step(struct engine *e) {
mask
|=
1
<<
task_type_grav_external
;
}
/* Add the tasks corresponding to cooling to the masks */
/* Add the tasks corresponding to cooling to the masks */
if
(
e
->
policy
&
engine_policy_cooling
)
{
mask
|=
1
<<
task_type_cooling
;
}
...
...
@@ -3115,8 +3113,7 @@ void engine_init(struct engine *e, struct space *s,
const
struct
phys_const
*
physical_constants
,
const
struct
hydro_props
*
hydro
,
const
struct
external_potential
*
potential
,
const
struct
cooling_data
*
cooling
)
{
const
struct
cooling_data
*
cooling
)
{
/* Clean-up everything */
bzero
(
e
,
sizeof
(
struct
engine
));
...
...
src/engine.h
View file @
47883b9f
...
...
@@ -222,7 +222,7 @@ void engine_init(struct engine *e, struct space *s,
const
struct
phys_const
*
physical_constants
,
const
struct
hydro_props
*
hydro
,
const
struct
external_potential
*
potential
,
const
struct
cooling_data
*
cooling
);
const
struct
cooling_data
*
cooling
);
void
engine_launch
(
struct
engine
*
e
,
int
nr_runners
,
unsigned
int
mask
,
unsigned
int
submask
);
void
engine_prepare
(
struct
engine
*
e
);
...
...
src/runner.c
View file @
47883b9f
...
...
@@ -147,7 +147,6 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
if
(
timer
)
TIMER_TOC
(
timer_dograv_external
);
}
/**
* @brief Calculate change in entropy from cooling
*
...
...
@@ -187,8 +186,8 @@ 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
(
constants
,
us
,
cooling
,
p
,
dt
);
dt
=
(
p
->
ti_end
-
p
->
ti_begin
)
*
timeBase
;
update_entropy
(
constants
,
us
,
cooling
,
p
,
dt
);
}
}
...
...
src/task.c
View file @
47883b9f
...
...
@@ -48,9 +48,9 @@
/* Task type names. */
const
char
*
taskID_names
[
task_type_count
]
=
{
"none"
,
"sort"
,
"self"
,
"pair"
,
"sub_self"
,
"sub_pair"
,
"init"
,
"ghost"
,
"extra_ghost"
,
"kick"
,
"kick_fixdt"
,
"send"
,
"recv"
,
"grav_gather_m"
,
"grav_fft"
,
"none"
,
"sort"
,
"self"
,
"pair"
,
"sub_self"
,
"sub_pair"
,
"init"
,
"ghost"
,
"extra_ghost"
,
"kick"
,
"kick_fixdt"
,
"send"
,
"recv"
,
"grav_gather_m"
,
"grav_fft"
,
"grav_mm"
,
"grav_up"
,
"grav_external"
,
"cooling"
};
const
char
*
subtaskID_names
[
task_subtype_count
]
=
{
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment