cooling.h 7.38 KB
Newer Older
Stefan Arridge's avatar
Stefan Arridge committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *                    Richard Bower (r.g.bower@durham.ac.uk)
 *                    Stefan Arridge  (stefan.arridge@durham.ac.uk)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

23
24
#ifndef SWIFT_COOLING_CONST_LAMBDA_H
#define SWIFT_COOLING_CONST_LAMBDA_H
Stefan Arridge's avatar
Stefan Arridge committed
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

/* Some standard headers. */
#include <math.h>

/* Local includes. */
#include "const.h"
#include "error.h"
#include "hydro.h"
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
#include "units.h"

/* Cooling Properties */
struct cooling_data {
40

Stefan Arridge's avatar
Stefan Arridge committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
  /*! Cooling rate in cgs units. Defined by 'rho * du/dt = -lambda * n_H^2'*/
  float lambda;

  /*! Minimum temperature (in Kelvin) for all gas particles*/
  float min_temperature;

  /*! Fraction of gas mass that is Hydrogen. Used to calculate n_H*/
  float hydrogen_mass_abundance;

  /* 'mu', used to convert min_temperature to min_internal energy*/
  float mean_molecular_weight;

  /*! Minimally allowed internal energy of the particles */
  float min_energy;
  float min_energy_cgs;

  /*! Constant multiplication factor for time-step criterion */
  float cooling_tstep_mult;
};

/**
 * @brief Calculates du/dt in code units for a particle.
 *
 * @param phys_const The physical constants in internal units.
 * @param us The internal system of units.
 * @param cooling The #cooling_data used in the run.
 * @param p Pointer to the particle data..
 */
__attribute__((always_inline)) INLINE static float cooling_rate(
    const struct phys_const* const phys_const, const struct UnitSystem* us,
    const struct cooling_data* cooling, const struct part* p) {

  /* Get particle properties */
  /* Density */
75
  const float rho = hydro_get_density(p);
Stefan Arridge's avatar
Stefan Arridge committed
76
77
78
79
80
81
  /* Get cooling function properties */
  const float X_H = cooling->hydrogen_mass_abundance;
  /* lambda should always be set in cgs units */
  const float lambda_cgs = cooling->lambda;

  /*convert from internal code units to cgs*/
82
83
  const float rho_cgs =
      rho * units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
Stefan Arridge's avatar
Stefan Arridge committed
84
  const float m_p_cgs = phys_const->const_proton_mass *
85
                        units_cgs_conversion_factor(us, UNIT_CONV_MASS);
Stefan Arridge's avatar
Stefan Arridge committed
86
87
88
89
90
91
  const float n_H_cgs = X_H * rho_cgs / m_p_cgs;

  /* Calculate du_dt */
  const float du_dt_cgs = -lambda_cgs * n_H_cgs * n_H_cgs / rho_cgs;

  /* Convert du/dt back to internal code units */
92
93
94
  const float du_dt =
      du_dt_cgs * units_cgs_conversion_factor(us, UNIT_CONV_TIME) /
      units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
Stefan Arridge's avatar
Stefan Arridge committed
95
96
97
98
99
100
101
102
103
104
105
106
107
108

  return du_dt;
}

/**
 * @brief Apply the cooling function to a particle.
 *
 * @param phys_const The physical constants in internal units.
 * @param us The internal system of units.
 * @param cooling The #cooling_data used in the run.
 * @param p Pointer to the particle data.
 * @param dt The time-step of this particle.
 */
__attribute__((always_inline)) INLINE static void cooling_cool_part(
109
110
111
112
    const struct phys_const* restrict phys_const,
    const struct UnitSystem* restrict us,
    const struct cooling_data* restrict cooling, struct part* restrict p,
    float dt) {
Stefan Arridge's avatar
Stefan Arridge committed
113
114
115
116
117
118
119
120

  /* 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;

  /* Calculate du_dt */
121
  const float du_dt = cooling_rate(phys_const, us, cooling, p);
122

Stefan Arridge's avatar
Stefan Arridge committed
123
124
125
126
127
128
129
  /* Intergrate cooling equation, but enforce energy floor */
  float u_new;
  if (u_old + du_dt * dt > u_floor) {
    u_new = u_old + du_dt * dt;
  } else {
    u_new = u_floor;
  }
130

Stefan Arridge's avatar
Stefan Arridge committed
131
132
  /* Update the internal energy */
  hydro_set_internal_energy(p, u_new);
133
134
135
136
137
138
139

  /* if (-(u_new_test - u_old) / u_old > 1.0e-6) */
  /*   error( */
  /*       "Particle has not successfully cooled: u_old = %g , du_dt = %g , dt =
   * " */
  /*       "%g, du_dt*dt = %g, u_old + du_dt*dt = %g, u_new = %g\n", */
  /*       u_old, du_dt, dt, du_dt * dt, u_new, u_new_test); */
Stefan Arridge's avatar
Stefan Arridge committed
140
141
142
143
144
145
146
147
}

/**
 * @brief Computes the time-step due to cooling
 *
 * @param cooling The #cooling_data used in the run.
 * @param phys_const The physical constants in internal units.
 * @param us The internal system of units.
Matthieu Schaller's avatar
Matthieu Schaller committed
148
 * @param p Pointer to the particle data.
Stefan Arridge's avatar
Stefan Arridge committed
149
150
 */
__attribute__((always_inline)) INLINE static float cooling_timestep(
151
152
153
    const struct cooling_data* restrict cooling,
    const struct phys_const* restrict phys_const,
    const struct UnitSystem* restrict us, const struct part* restrict p) {
Stefan Arridge's avatar
Stefan Arridge committed
154
155

  /* Get du_dt */
156
157
  const float du_dt = cooling_rate(phys_const, us, cooling, p);

Stefan Arridge's avatar
Stefan Arridge committed
158
159
160
  /* Get current internal energy (dt=0) */
  const float u = hydro_get_internal_energy(p, 0.f);

161
  return u / du_dt;
Stefan Arridge's avatar
Stefan Arridge committed
162
163
164
165
166
167
168
169
170
171
}

/**
 * @brief Initialises the cooling properties.
 *
 * @param parameter_file The parsed parameter file.
 * @param us The current internal system of units.
 * @param phys_const The physical constants in internal units.
 * @param cooling The cooling properties to initialize
 */
172
173
174
static INLINE void cooling_init_backend(
    const struct swift_params* parameter_file, const struct UnitSystem* us,
    const struct phys_const* phys_const, struct cooling_data* cooling) {
Stefan Arridge's avatar
Stefan Arridge committed
175

176
177
178
179
  cooling->lambda =
      parser_get_param_double(parameter_file, "LambdaCooling:lambda");
  cooling->min_temperature = parser_get_param_double(
      parameter_file, "LambdaCooling:minimum_temperature");
180
  cooling->hydrogen_mass_abundance = parser_get_param_double(
181
182
183
184
185
      parameter_file, "LambdaCooling:hydrogen_mass_abundance");
  cooling->mean_molecular_weight = parser_get_param_double(
      parameter_file, "LambdaCooling:mean_molecular_weight");
  cooling->cooling_tstep_mult = parser_get_param_double(
      parameter_file, "LambdaCooling:cooling_tstep_mult");
Stefan Arridge's avatar
Stefan Arridge committed
186
187
188

  /*convert minimum temperature into minimum internal energy*/
  const float u_floor =
189
190
191
      phys_const->const_boltzmann_k * cooling->min_temperature /
      (hydro_gamma_minus_one * cooling->mean_molecular_weight *
       phys_const->const_proton_mass);
Stefan Arridge's avatar
Stefan Arridge committed
192
193
194
195
196
197
198
199
200
201
202
203
  const float u_floor_cgs =
      u_floor * units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);

  cooling->min_energy = u_floor;
  cooling->min_energy_cgs = u_floor_cgs;
}

/**
 * @brief Prints the properties of the cooling model to stdout.
 *
 * @param cooling The properties of the cooling function.
 */
204
static INLINE void cooling_print_backend(const struct cooling_data* cooling) {
Stefan Arridge's avatar
Stefan Arridge committed
205

206
207
208
209
210
211
  message(
      "Cooling function is 'Constant lambda' with "
      "(lambda,min_temperature,hydrogen_mass_abundance,mean_molecular_weight) "
      "=  (%g,%g,%g,%g)",
      cooling->lambda, cooling->min_temperature,
      cooling->hydrogen_mass_abundance, cooling->mean_molecular_weight);
Stefan Arridge's avatar
Stefan Arridge committed
212
213
}

214
#endif /* SWIFT_COOLING_CONST_LAMBDA_H */