cooling.h 12 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@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/>.
 *
 ******************************************************************************/
#ifndef SWIFT_COOLING_GRACKLE_H
#define SWIFT_COOLING_GRACKLE_H

/**
 * @file src/cooling/none/cooling.h
 * @brief Empty infrastructure for the cases without cooling function
 */

/* Some standard headers. */
#include <float.h>
29
#include <grackle.h>
Matthieu Schaller's avatar
Matthieu Schaller committed
30
#include <math.h>
31
32
33
34

/* Local includes. */
#include "error.h"
#include "hydro.h"
lhausamm's avatar
lhausamm committed
35
#include "io_properties.h"
36
37
38
39
40
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
#include "units.h"

lhausamm's avatar
lhausamm committed
41
/* need to rework (and check) code if changed */
42
#define GRACKLE_NPART 1
lhausamm's avatar
lhausamm committed
43
#define GRACKLE_RANK 3
44

45
46
#ifdef HAVE_HDF5

lhausamm's avatar
lhausamm committed
47
48
49
50
/**
 * @brief Writes the current model of SPH to the file
 * @param h_grpsph The HDF5 group in which to write
 */
51
__attribute__((always_inline)) INLINE static void cooling_write_flavour(
52
    hid_t h_grpsph) {
lhausamm's avatar
lhausamm committed
53

54
  io_write_attribute_s(h_grpsph, "Cooling Model", "Grackle");
lhausamm's avatar
lhausamm committed
55
}
56
#endif
lhausamm's avatar
lhausamm committed
57

58
/**
59
60
 * @brief Sets the cooling properties of the (x-)particles to a valid start
 * state.
61
 *
62
63
 * @param p Pointer to the particle data.
 * @param xp Pointer to the extended particle data.
64
 * @param cooling The properties of the cooling function.
65
 */
66
67
68
__attribute__((always_inline)) INLINE static void cooling_first_init_part(
    const struct part* restrict p, struct xpart* restrict xp,
    const struct cooling_function_data* cooling) {
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83

  xp->cooling_data.radiated_energy = 0.f;
}

/**
 * @brief Returns the total radiated energy by this particle.
 *
 * @param xp The extended particle data
 */
__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
    const struct xpart* restrict xp) {

  return xp->cooling_data.radiated_energy;
}

84
85
86
87
88
/**
 * @brief Prints the properties of the cooling model to stdout.
 *
 * @param cooling The properties of the cooling function.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
89
__attribute__((always_inline)) INLINE static void cooling_print_backend(
90
91
92
    const struct cooling_function_data* cooling) {

  message("Cooling function is 'Grackle'.");
93
  message("Using Grackle           = %i", cooling->chemistry.use_grackle);
Matthieu Schaller's avatar
Matthieu Schaller committed
94
95
96
97
  message("Chemical network        = %i",
          cooling->chemistry.primordial_chemistry);
  message("Radiative cooling       = %i",
          cooling->chemistry.with_radiative_cooling);
98
  message("Metal cooling           = %i", cooling->chemistry.metal_cooling);
Matthieu Schaller's avatar
Matthieu Schaller committed
99
100

  message("CloudyTable             = %s", cooling->cloudy_table);
101
102
  message("UVbackground            = %d", cooling->uv_background);
  message("Redshift                = %g", cooling->redshift);
Matthieu Schaller's avatar
Matthieu Schaller committed
103
  message("Density Self Shielding  = %g", cooling->density_self_shielding);
104
105
  message("Units:");
  message("\tComoving     = %i", cooling->units.comoving_coordinates);
lhausamm's avatar
lhausamm committed
106
107
  message("\tLength       = %g", cooling->units.length_units);
  message("\tDensity      = %g", cooling->units.density_units);
108
  message("\tTime         = %g", cooling->units.time_units);
lhausamm's avatar
lhausamm committed
109
  message("\tScale Factor = %g", cooling->units.a_units);
110
111
}

112
113
/**
 * @brief Compute the cooling rate
114
115
116
117
118
119
 *
 * @param phys_const The physical constants in internal units.
 * @param us The internal system of units.
 * @param cooling The #cooling_function_data used in the run.
 * @param p Pointer to the particle data.
 * @param dt The time-step of this particle.
120
121
 *
 * @return du / dt
122
 */
123
124
125
__attribute__((always_inline)) INLINE static double cooling_rate(
    const struct phys_const* restrict phys_const,
    const struct unit_system* restrict us,
126
    const struct cosmology* restrict cosmo,
127
128
    const struct cooling_function_data* restrict cooling,
    struct part* restrict p, float dt) {
129

Matthieu Schaller's avatar
Matthieu Schaller committed
130
  if (cooling->chemistry.primordial_chemistry > 1) error("Not implemented");
lhausamm's avatar
lhausamm committed
131

lhausamm's avatar
lhausamm committed
132
  /* set current time */
lhausamm's avatar
lhausamm committed
133
  code_units units = cooling->units;
lhausamm's avatar
lhausamm committed
134
135
136
  if (cooling->redshift == -1)
    error("TODO time dependant redshift");
  else
lhausamm's avatar
lhausamm committed
137
    units.a_value = 1. / (1. + cooling->redshift);
138

lhausamm's avatar
lhausamm committed
139
140
  /* initialize data */
  grackle_field_data data;
lhausamm's avatar
lhausamm committed
141

lhausamm's avatar
lhausamm committed
142
143
144
145
146
  /* set values */
  /* grid */
  int grid_dimension[GRACKLE_RANK] = {GRACKLE_NPART, 1, 1};
  int grid_start[GRACKLE_RANK] = {0, 0, 0};
  int grid_end[GRACKLE_RANK] = {GRACKLE_NPART - 1, 0, 0};
Matthieu Schaller's avatar
Matthieu Schaller committed
147

lhausamm's avatar
lhausamm committed
148
149
150
151
  data.grid_rank = GRACKLE_RANK;
  data.grid_dimension = grid_dimension;
  data.grid_start = grid_start;
  data.grid_end = grid_end;
lhausamm's avatar
lhausamm committed
152

lhausamm's avatar
lhausamm committed
153
  /* general particle data */
154
155
  gr_float density = hydro_get_physical_density(p, cosmo);
  const double energy_before = hydro_get_physical_internal_energy(p, cosmo);
lhausamm's avatar
lhausamm committed
156
157
158
159
160
161
162
163
164
165
166
  gr_float energy = energy_before;
  gr_float vx = 0;
  gr_float vy = 0;
  gr_float vz = 0;

  data.density = &density;
  data.internal_energy = &energy;
  data.x_velocity = &vx;
  data.y_velocity = &vy;
  data.z_velocity = &vz;

167
168
169
170
171
172
173
  /* /\* primordial chemistry >= 1 *\/ */
  /* gr_float HI_density = density; */
  /* gr_float HII_density = 0.; */
  /* gr_float HeI_density = 0.; */
  /* gr_float HeII_density = 0.; */
  /* gr_float HeIII_density = 0.; */
  /* gr_float e_density = 0.; */
Matthieu Schaller's avatar
Matthieu Schaller committed
174

175
176
177
178
179
180
181
182
183
184
185
  /* data.HI_density = &HI_density; */
  /* data.HII_density = &HII_density; */
  /* data.HeI_density = &HeI_density; */
  /* data.HeII_density = &HeII_density; */
  /* data.HeIII_density = &HeIII_density; */
  /* data.e_density = &e_density; */

  /* /\* primordial chemistry >= 2 *\/ */
  /* gr_float HM_density = 0.; */
  /* gr_float H2I_density = 0.; */
  /* gr_float H2II_density = 0.; */
Matthieu Schaller's avatar
Matthieu Schaller committed
186

187
188
189
  /* data.HM_density = &HM_density; */
  /* data.H2I_density = &H2I_density; */
  /* data.H2II_density = &H2II_density; */
Matthieu Schaller's avatar
Matthieu Schaller committed
190

191
192
193
194
  /* /\* primordial chemistry >= 3 *\/ */
  /* gr_float DI_density = 0.; */
  /* gr_float DII_density = 0.; */
  /* gr_float HDI_density = 0.; */
Matthieu Schaller's avatar
Matthieu Schaller committed
195

196
197
198
199
  /* data.DI_density = &DI_density; */
  /* data.DII_density = &DII_density; */
  /* data.HDI_density = &HDI_density; */

lhausamm's avatar
lhausamm committed
200
  /* metal cooling = 1 */
201
  gr_float metal_density = density * grackle_data->SolarMetalFractionByMass;
lhausamm's avatar
lhausamm committed
202

lhausamm's avatar
lhausamm committed
203
  data.metal_density = &metal_density;
lhausamm's avatar
lhausamm committed
204

205
206
207
208
  /* /\* volumetric heating rate *\/ */
  /* gr_float volumetric_heating_rate = 0.; */

  /* data.volumetric_heating_rate = &volumetric_heating_rate; */
Matthieu Schaller's avatar
Matthieu Schaller committed
209

210
211
212
213
214
  /* /\* specific heating rate *\/ */
  /* gr_float specific_heating_rate = 0.; */

  /* data.specific_heating_rate = &specific_heating_rate; */

lhausamm's avatar
lhausamm committed
215
  /* solve chemistry with table */
lhausamm's avatar
lhausamm committed
216
  if (solve_chemistry(&units, &data, dt) == 0) {
lhausamm's avatar
lhausamm committed
217
    error("Error in solve_chemistry.");
218
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
219

lhausamm's avatar
lhausamm committed
220
  return (energy - energy_before) / dt;
221
222
223
224
225
226
227
}

/**
 * @brief Apply the cooling function to a particle.
 *
 * @param phys_const The physical constants in internal units.
 * @param us The internal system of units.
228
 * @param cosmo The current cosmological model.
229
230
231
232
233
234
 * @param cooling The #cooling_function_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(
    const struct phys_const* restrict phys_const,
lhausamm's avatar
lhausamm committed
235
    const struct unit_system* restrict us,
236
    const struct cosmology* restrict cosmo,
237
238
239
    const struct cooling_function_data* restrict cooling,
    struct part* restrict p, struct xpart* restrict xp, float dt) {

Matthieu Schaller's avatar
Matthieu Schaller committed
240
  if (dt == 0.) return;
241
242
243

  /* Current du_dt */
  const float hydro_du_dt = hydro_get_internal_energy_dt(p);
244

245
  /* compute cooling rate */
246
  const float du_dt = cooling_rate(phys_const, us, cosmo, cooling, p, dt);
247
248

  /* record energy lost */
Matthieu Schaller's avatar
Matthieu Schaller committed
249
  xp->cooling_data.radiated_energy += -du_dt * dt * hydro_get_mass(p);
lhausamm's avatar
lhausamm committed
250

251
  /* Update the internal energy */
252
  hydro_set_internal_energy_dt(p, hydro_du_dt + du_dt);
253
254
255
256
257
258
259
260
261
}

/**
 * @brief Computes the cooling time-step.
 *
 * We return FLT_MAX so as to impose no limit on the time-step.
 *
 * @param cooling The #cooling_function_data used in the run.
 * @param phys_const The physical constants in internal units.
262
 * @param cosmo The current cosmological model.
263
264
265
266
267
268
 * @param us The internal system of units.
 * @param p Pointer to the particle data.
 */
__attribute__((always_inline)) INLINE static float cooling_timestep(
    const struct cooling_function_data* restrict cooling,
    const struct phys_const* restrict phys_const,
269
    const struct cosmology* restrict cosmo,
lhausamm's avatar
lhausamm committed
270
    const struct unit_system* restrict us, const struct part* restrict p) {
271
272
273
274
275
276
277
278
279
280
281
282

  return FLT_MAX;
}

/**
 * @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
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
283
__attribute__((always_inline)) INLINE static void cooling_init_backend(
lhausamm's avatar
lhausamm committed
284
    const struct swift_params* parameter_file, const struct unit_system* us,
285
286
287
    const struct phys_const* phys_const,
    struct cooling_function_data* cooling) {

Matthieu Schaller's avatar
Matthieu Schaller committed
288
  /* read parameters */
lhausamm's avatar
lhausamm committed
289
  parser_get_param_string(parameter_file, "GrackleCooling:GrackleCloudyTable",
290
291
                          cooling->cloudy_table);
  cooling->uv_background =
Matthieu Schaller's avatar
Matthieu Schaller committed
292
      parser_get_param_int(parameter_file, "GrackleCooling:UVbackground");
293
294

  cooling->redshift =
Matthieu Schaller's avatar
Matthieu Schaller committed
295
      parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift");
296
297

  cooling->density_self_shielding = parser_get_param_double(
lhausamm's avatar
lhausamm committed
298
      parameter_file, "GrackleCooling:GrackleHSShieldingDensityThreshold");
Matthieu Schaller's avatar
Matthieu Schaller committed
299

300
301
302
303
#ifdef SWIFT_DEBUG_CHECKS
  /* enable verbose for grackle */
  grackle_verbose = 1;
#endif
lhausamm's avatar
lhausamm committed
304

305
306
307
308
309
  /* Set up the units system.
     These are conversions from code units to cgs. */

  /* first cosmo */
  cooling->units.a_units = 1.0;  // units for the expansion factor (1/1+zi)
310
  cooling->units.a_value = 1.0;
311
312
313

  /* We assume here all physical quantities to
     be in proper coordinate (not comobile)  */
314
  cooling->units.comoving_coordinates = 0;
315
316

  /* then units */
Matthieu Schaller's avatar
Matthieu Schaller committed
317
318
  cooling->units.density_units =
      us->UnitMass_in_cgs / pow(us->UnitLength_in_cgs, 3);
319
320
  cooling->units.length_units = us->UnitLength_in_cgs;
  cooling->units.time_units = us->UnitTime_in_cgs;
Matthieu Schaller's avatar
Matthieu Schaller committed
321
322
323
324
325
  cooling->units.velocity_units = cooling->units.a_units *
                                  cooling->units.length_units /
                                  cooling->units.time_units;

  chemistry_data* chemistry = &cooling->chemistry;
lhausamm's avatar
lhausamm committed
326

327
  /* Create a chemistry object for parameters and rate data. */
lhausamm's avatar
lhausamm committed
328
  if (set_default_chemistry_parameters(chemistry) == 0) {
329
330
    error("Error in set_default_chemistry_parameters.");
  }
lhausamm's avatar
lhausamm committed
331

332
  // Set parameter values for chemistry.
lhausamm's avatar
lhausamm committed
333
334
  chemistry->use_grackle = 1;
  chemistry->with_radiative_cooling = 1;
335
336
  /* molecular network with H, He, D
   From Cloudy table */
lhausamm's avatar
lhausamm committed
337
338
339
340
  chemistry->primordial_chemistry = 0;
  chemistry->metal_cooling = 1;  // metal cooling on
  chemistry->UVbackground = cooling->uv_background;
  chemistry->grackle_data_file = cooling->cloudy_table;
341
342
343
  chemistry->use_radiative_transfer = 0;
  chemistry->use_volumetric_heating_rate = 0;
  chemistry->use_specific_heating_rate = 0;
344

lhausamm's avatar
lhausamm committed
345
346
  /* Initialize the chemistry object. */
  if (initialize_chemistry_data(&cooling->units) == 0) {
347
348
    error("Error in initialize_chemistry_data.");
  }
lhausamm's avatar
lhausamm committed
349
350

#ifdef SWIFT_DEBUG_CHECKS
351
352
  if (GRACKLE_NPART != 1)
    error("Grackle with multiple particles not implemented");
353
  float threshold = cooling->density_self_shielding;
lhausamm's avatar
lhausamm committed
354
355
356
357
358
359
360

  threshold /= phys_const->const_proton_mass;
  threshold /= pow(us->UnitLength_in_cgs, 3);

  message("***************************************");
  message("initializing grackle cooling function");
  message("");
361
362
  cooling_print_backend(cooling);
  message("Density Self Shielding = %g atom/cm3", threshold);
lhausamm's avatar
lhausamm committed
363
364
365
366

  message("");
  message("***************************************");
#endif
367
368
369
}

#endif /* SWIFT_COOLING_GRACKLE_H */