cooling.h 10.4 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'.");
lhausamm's avatar
lhausamm committed
78
79
80
81
82
  message("Using Grackle           = %i", grackle_data->use_grackle);
  message("Chemical network        = %i", grackle_data->primordial_chemistry);
  message("Radiative cooling       = %i", grackle_data->with_radiative_cooling);
  message("Metal cooling           = %i", grackle_data->metal_cooling);
  
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

lhausamm's avatar
lhausamm committed
116
117
  cooling_print_backend(cooling);

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

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

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

lhausamm's avatar
lhausamm committed
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
  /* 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;

  /* primordial chemistry >= 1 */
  data.HI_density = NULL;
  data.HII_density = NULL;
  data.HeI_density = NULL;
  data.HeII_density = NULL;
  data.HeIII_density = NULL;
  data.e_density = NULL;
  /* primordial chemistry >= 2 */
  data.HM_density = NULL;
  data.H2I_density = NULL;
  data.H2II_density = NULL;
  /* primordial chemistry >= 3 */
  data.DI_density = NULL;
  data.DII_density = NULL;
  data.HDI_density = NULL;
  /* metal cooling = 1 */
  const double Z = 0.02041;
  gr_float metal_density = Z * density;
lhausamm's avatar
lhausamm committed
171

lhausamm's avatar
lhausamm committed
172
  data.metal_density = &metal_density;
lhausamm's avatar
lhausamm committed
173

lhausamm's avatar
lhausamm committed
174
175
176
177
178
  /* volumetric heating rate */
  data.volumetric_heating_rate = NULL;
  /* specific heating rate */
  data.specific_heating_rate = NULL;
  
lhausamm's avatar
lhausamm committed
179
  /* solve chemistry with table */
lhausamm's avatar
lhausamm committed
180
  if (solve_chemistry(&units, &data, dt) == 0) {
lhausamm's avatar
lhausamm committed
181
    error("Error in solve_chemistry.");
182
  }
lhausamm's avatar
lhausamm committed
183
  message("Energy: %g, %g, %g", energy, energy_before, dt);
184

lhausamm's avatar
lhausamm committed
185
186
  exit(1);
  return (energy - energy_before) / dt;
187
188
189
190
191
192
193
194
195
196
197
198
199
}

/**
 * @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
200
    const struct unit_system* restrict us,
201
202
203
    const struct cooling_function_data* restrict cooling,
    struct part* restrict p, struct xpart* restrict xp, float dt) {

Matthieu Schaller's avatar
Matthieu Schaller committed
204
  if (dt == 0.) return;
205
206
207

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

209
210
  /* compute cooling rate */
  const float du_dt = cooling_rate(phys_const, us, cooling, p, dt);
211
212

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

215
  /* Update the internal energy */
216
  hydro_set_internal_energy_dt(p, hydro_du_dt + du_dt);
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
}

/**
 * @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
232
    const struct unit_system* restrict us, const struct part* restrict p) {
233
234
235
236
237
238
239
240
241
242
243
244

  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
 */
245
__attribute__((always_inline))INLINE static void cooling_init_backend(
lhausamm's avatar
lhausamm committed
246
    const struct swift_params* parameter_file, const struct unit_system* us,
247
248
249
    const struct phys_const* phys_const,
    struct cooling_function_data* cooling) {

250
    /* read parameters */
lhausamm's avatar
lhausamm committed
251
  parser_get_param_string(parameter_file, "GrackleCooling:GrackleCloudyTable",
252
253
254
255
256
257
258
259
                          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
260
      parameter_file, "GrackleCooling:GrackleHSShieldingDensityThreshold");
261
262
263
264
265
  
#ifdef SWIFT_DEBUG_CHECKS
  /* enable verbose for grackle */
  grackle_verbose = 1;
#endif
lhausamm's avatar
lhausamm committed
266

267
268
269
270
271
272
273
274
  /* 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)

  /* We assume here all physical quantities to
     be in proper coordinate (not comobile)  */
275
  cooling->units.comoving_coordinates = 0;
276
277
278
279
280
281
282

  /* 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
283
284
285

  chemistry_data *chemistry = &cooling->chemistry;
  
286
  /* Create a chemistry object for parameters and rate data. */
lhausamm's avatar
lhausamm committed
287
  if (set_default_chemistry_parameters(chemistry) == 0) {
288
289
    error("Error in set_default_chemistry_parameters.");
  }
lhausamm's avatar
lhausamm committed
290

291
  // Set parameter values for chemistry.
lhausamm's avatar
lhausamm committed
292
293
  chemistry->use_grackle = 1;
  chemistry->with_radiative_cooling = 1;
294
295
  /* molecular network with H, He, D
   From Cloudy table */
lhausamm's avatar
lhausamm committed
296
297
298
299
  chemistry->primordial_chemistry = 0;
  chemistry->metal_cooling = 1;  // metal cooling on
  chemistry->UVbackground = cooling->uv_background;
  chemistry->grackle_data_file = cooling->cloudy_table;
300

lhausamm's avatar
lhausamm committed
301
302
  /* Initialize the chemistry object. */
  if (initialize_chemistry_data(&cooling->units) == 0) {
303
304
    error("Error in initialize_chemistry_data.");
  }
lhausamm's avatar
lhausamm committed
305
306

#ifdef SWIFT_DEBUG_CHECKS
307
308
  if (GRACKLE_NPART != 1)
    error("Grackle with multiple particles not implemented");
lhausamm's avatar
lhausamm committed
309
310
311
312
313
314
315
316
  float threshold = cooling->GrackleHSShieldingDensityThreshold;

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

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

320

321
  //grackle_print_data();
lhausamm's avatar
lhausamm committed
322
323
324
  message("");
  message("***************************************");
#endif
325
326
327
}

#endif /* SWIFT_COOLING_GRACKLE_H */