cooling.h 11.3 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
29
/*******************************************************************************
 * 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>
#include <math.h>
30
#include <grackle.h>
31
32
33
34
35
36
37
38
39

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

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

44
/**
45
46
 * @brief Sets the cooling properties of the (x-)particles to a valid start
 * state.
47
 *
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
 * @param p Pointer to the particle data.
 * @param xp Pointer to the extended particle data.
 */
__attribute__((always_inline)) INLINE static void cooling_init_part(
    const struct part* restrict p, struct xpart* restrict xp) {

  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;
}


69
70
71
72
73
74
75
76
77
/**
 * @brief Prints the properties of the cooling model to stdout.
 *
 * @param cooling The properties of the cooling function.
 */
__attribute__((always_inline))INLINE static void cooling_print_backend(
    const struct cooling_function_data* cooling) {

  message("Cooling function is 'Grackle'.");
78
79
80
81
  message("Using Grackle           = %i", cooling->chemistry.use_grackle);
  message("Chemical network        = %i", cooling->chemistry.primordial_chemistry);
  message("Radiative cooling       = %i", cooling->chemistry.with_radiative_cooling);
  message("Metal cooling           = %i", cooling->chemistry.metal_cooling);
lhausamm's avatar
lhausamm committed
82
  
83
84
85
86
87
88
89
90
  message("CloudyTable             = %s",
          cooling->cloudy_table);
  message("UVbackground            = %d", cooling->uv_background);
  message("Redshift                = %g", cooling->redshift);
  message("Density Self Shielding  = %g",
          cooling->density_self_shielding);
  message("Units:");
  message("\tComoving     = %i", cooling->units.comoving_coordinates);
lhausamm's avatar
lhausamm committed
91
92
  message("\tLength       = %g", cooling->units.length_units);
  message("\tDensity      = %g", cooling->units.density_units);
93
  message("\tTime         = %g", cooling->units.time_units);
lhausamm's avatar
lhausamm committed
94
95
  message("\tScale Factor = %g", cooling->units.a_units);

96
97
}

lhausamm's avatar
lhausamm committed
98

99
100
/**
 * @brief Compute the cooling rate
101
102
103
104
105
106
 *
 * @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.
107
108
 *
 * @return du / dt
109
 */
110
111
112
113
114
__attribute__((always_inline)) INLINE static double cooling_rate(
    const struct phys_const* restrict phys_const,
    const struct unit_system* restrict us,
    const struct cooling_function_data* restrict cooling,
    struct part* restrict p, float dt) {
115

116
117
  if (cooling->chemistry.primordial_chemistry > 1)
    error("Not implemented");
lhausamm's avatar
lhausamm committed
118

lhausamm's avatar
lhausamm committed
119
  /* set current time */
lhausamm's avatar
lhausamm committed
120
  code_units units = cooling->units;
lhausamm's avatar
lhausamm committed
121
122
123
  if (cooling->redshift == -1)
    error("TODO time dependant redshift");
  else
lhausamm's avatar
lhausamm committed
124
    units.a_value = 1. / (1. + cooling->redshift);
125

lhausamm's avatar
lhausamm committed
126
127
  /* initialize data */
  grackle_field_data data;
lhausamm's avatar
lhausamm committed
128

lhausamm's avatar
lhausamm committed
129
130
131
132
133
134
135
136
137
138
  /* 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};
  
  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
139

lhausamm's avatar
lhausamm committed
140
141
142
143
144
145
146
147
148
149
150
151
152
153
  /* general particle data */
  gr_float density = hydro_get_density(p);
  const double energy_before = hydro_get_internal_energy(p);
  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;

154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
  /* /\* 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.; */
  
  /* 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.; */
  
  /* data.HM_density = &HM_density; */
  /* data.H2I_density = &H2I_density; */
  /* data.H2II_density = &H2II_density; */
  
  /* /\* primordial chemistry >= 3 *\/ */
  /* gr_float DI_density = 0.; */
  /* gr_float DII_density = 0.; */
  /* gr_float HDI_density = 0.; */
  
  /* data.DI_density = &DI_density; */
  /* data.DII_density = &DII_density; */
  /* data.HDI_density = &HDI_density; */

lhausamm's avatar
lhausamm committed
187
  /* metal cooling = 1 */
188
  gr_float metal_density = density * grackle_data->SolarMetalFractionByMass;
lhausamm's avatar
lhausamm committed
189

lhausamm's avatar
lhausamm committed
190
  data.metal_density = &metal_density;
lhausamm's avatar
lhausamm committed
191

192
193
194
195
  /* /\* volumetric heating rate *\/ */
  /* gr_float volumetric_heating_rate = 0.; */

  /* data.volumetric_heating_rate = &volumetric_heating_rate; */
lhausamm's avatar
lhausamm committed
196
  
197
198
199
200
201
  /* /\* specific heating rate *\/ */
  /* gr_float specific_heating_rate = 0.; */

  /* data.specific_heating_rate = &specific_heating_rate; */

lhausamm's avatar
lhausamm committed
202
  /* solve chemistry with table */
lhausamm's avatar
lhausamm committed
203
  if (solve_chemistry(&units, &data, dt) == 0) {
lhausamm's avatar
lhausamm committed
204
    error("Error in solve_chemistry.");
205
  }
206
  
lhausamm's avatar
lhausamm committed
207
  return (energy - energy_before) / dt;
208
209
210
211
212
213
214
215
216
217
218
219
220
}

/**
 * @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_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
221
    const struct unit_system* restrict us,
222
223
224
    const struct cooling_function_data* restrict cooling,
    struct part* restrict p, struct xpart* restrict xp, float dt) {

Matthieu Schaller's avatar
Matthieu Schaller committed
225
  if (dt == 0.) return;
226
227
228

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

230
231
  /* compute cooling rate */
  const float du_dt = cooling_rate(phys_const, us, cooling, p, dt);
232
233

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

236
  /* Update the internal energy */
237
  hydro_set_internal_energy_dt(p, hydro_du_dt + du_dt);
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
}

/**
 * @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.
 * @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,
lhausamm's avatar
lhausamm committed
253
    const struct unit_system* restrict us, const struct part* restrict p) {
254
255
256
257
258
259
260
261
262
263
264
265

  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
 */
266
__attribute__((always_inline))INLINE static void cooling_init_backend(
lhausamm's avatar
lhausamm committed
267
    const struct swift_params* parameter_file, const struct unit_system* us,
268
269
270
    const struct phys_const* phys_const,
    struct cooling_function_data* cooling) {

271
    /* read parameters */
lhausamm's avatar
lhausamm committed
272
  parser_get_param_string(parameter_file, "GrackleCooling:GrackleCloudyTable",
273
274
275
276
277
278
279
280
                          cooling->cloudy_table);
  cooling->uv_background =
    parser_get_param_int(parameter_file, "GrackleCooling:UVbackground");

  cooling->redshift =
    parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift");

  cooling->density_self_shielding = parser_get_param_double(
lhausamm's avatar
lhausamm committed
281
      parameter_file, "GrackleCooling:GrackleHSShieldingDensityThreshold");
282
283
284
285
286
  
#ifdef SWIFT_DEBUG_CHECKS
  /* enable verbose for grackle */
  grackle_verbose = 1;
#endif
lhausamm's avatar
lhausamm committed
287

288
289
290
291
292
  /* 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)
293
  cooling->units.a_value = 1.0;
294
295
296

  /* We assume here all physical quantities to
     be in proper coordinate (not comobile)  */
297
  cooling->units.comoving_coordinates = 0;
298
299
300
301
302
303
304

  /* then units */
  cooling->units.density_units = us->UnitMass_in_cgs / pow(us->UnitLength_in_cgs, 3);
  cooling->units.length_units = us->UnitLength_in_cgs;
  cooling->units.time_units = us->UnitTime_in_cgs;
  cooling->units.velocity_units =
    cooling->units.a_units * cooling->units.length_units / cooling->units.time_units;
lhausamm's avatar
lhausamm committed
305
306
307

  chemistry_data *chemistry = &cooling->chemistry;
  
308
  /* Create a chemistry object for parameters and rate data. */
lhausamm's avatar
lhausamm committed
309
  if (set_default_chemistry_parameters(chemistry) == 0) {
310
311
    error("Error in set_default_chemistry_parameters.");
  }
lhausamm's avatar
lhausamm committed
312

313
  // Set parameter values for chemistry.
lhausamm's avatar
lhausamm committed
314
315
  chemistry->use_grackle = 1;
  chemistry->with_radiative_cooling = 1;
316
317
  /* molecular network with H, He, D
   From Cloudy table */
lhausamm's avatar
lhausamm committed
318
319
320
321
  chemistry->primordial_chemistry = 0;
  chemistry->metal_cooling = 1;  // metal cooling on
  chemistry->UVbackground = cooling->uv_background;
  chemistry->grackle_data_file = cooling->cloudy_table;
322
323
324
  chemistry->use_radiative_transfer = 0;
  chemistry->use_volumetric_heating_rate = 0;
  chemistry->use_specific_heating_rate = 0;
325

lhausamm's avatar
lhausamm committed
326
327
  /* Initialize the chemistry object. */
  if (initialize_chemistry_data(&cooling->units) == 0) {
328
329
    error("Error in initialize_chemistry_data.");
  }
lhausamm's avatar
lhausamm committed
330

331

lhausamm's avatar
lhausamm committed
332
#ifdef SWIFT_DEBUG_CHECKS
333
334
  if (GRACKLE_NPART != 1)
    error("Grackle with multiple particles not implemented");
335
  float threshold = cooling->density_self_shielding;
lhausamm's avatar
lhausamm committed
336
337
338
339
340
341
342

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

  message("***************************************");
  message("initializing grackle cooling function");
  message("");
343
344
  cooling_print_backend(cooling);
  message("Density Self Shielding = %g atom/cm3", threshold);
lhausamm's avatar
lhausamm committed
345

346

lhausamm's avatar
lhausamm committed
347
348
349
  message("");
  message("***************************************");
#endif
350

351
352
353
}

#endif /* SWIFT_COOLING_GRACKLE_H */