cooling.h 11.7 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
42
#define GRACKLE_NPART 1

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
69
70
 * @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;
}


/**
 * @brief Compute the cooling rate
71
72
73
74
75
76
 *
 * @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.
77
78
 *
 * @return du / dt
79
 */
80
81
82
83
84
__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) {
85

lhausamm's avatar
lhausamm committed
86
87
88
89
90
91
  /* set current time */
  float scale_factor;
  if (cooling->redshift == -1)
    error("TODO time dependant redshift");
  else
    scale_factor = 1. / (1. + cooling->redshift);
92

93
  /* Get current internal energy (dt=0) */
lhausamm's avatar
lhausamm committed
94
95
  const double energy_before = hydro_get_internal_energy(p);
  
96
97
  /* Get current density */
  const float rho = hydro_get_density(p);
lhausamm's avatar
lhausamm committed
98
  
99
100
  /* 0.02041 (= 1 Zsun in Grackle v2.0, but = 1.5761 Zsun in
     Grackle v2.1) */
lhausamm's avatar
lhausamm committed
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
  const double Z = 0.02041;

  /* create grackle struct */
  /* velocities */
  gr_float x_velocity[GRACKLE_NPART] = {0.0};
  gr_float y_velocity[GRACKLE_NPART] = {0.0};
  gr_float z_velocity[GRACKLE_NPART] = {0.0};

  /* particle data */
  gr_float density[GRACKLE_NPART] = {rho};
  gr_float metal_density[GRACKLE_NPART] = {Z * density[0]};
  gr_float energy[GRACKLE_NPART] = {energy_before};

  /* dimensions */
  int grid_dimension[3] = {GRACKLE_NPART, 0, 0};
  int grid_start[3] = {0, 0, 0};
  int grid_end[3] = {0, 0, 0};

  /* solve chemistry with table */
  if (solve_chemistry_table(&cooling->units, scale_factor, dt, grid_rank, grid_dimension,
                            grid_start, grid_end, density, energy, x_velocity,
                            y_velocity, z_velocity, metal_density) == 0) {
    error("Error in solve_chemistry.");
124
125
  }

lhausamm's avatar
lhausamm committed
126
  return (energy[0] - energy_before) / dt;
127
128
129
130
131
132
133
134
135
136
137
138
139
}

/**
 * @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
140
    const struct unit_system* restrict us,
141
142
143
    const struct cooling_function_data* restrict cooling,
    struct part* restrict p, struct xpart* restrict xp, float dt) {

Matthieu Schaller's avatar
Matthieu Schaller committed
144
  if (dt == 0.) return;
145
146
147

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

149
150
  /* compute cooling rate */
  const float du_dt = cooling_rate(phys_const, us, cooling, p, dt);
151
152

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

155
  /* Update the internal energy */
156
  hydro_set_internal_energy_dt(p, hydro_du_dt + du_dt);
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
}

/**
 * @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
172
    const struct unit_system* restrict us, const struct part* restrict p) {
173
174
175
176
177
178
179
180
181
182
183
184
185

  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
 */
static INLINE void cooling_init_backend(
lhausamm's avatar
lhausamm committed
186
    const struct swift_params* parameter_file, const struct unit_system* us,
187
188
189
    const struct phys_const* phys_const,
    struct cooling_function_data* cooling) {

190
    /* read parameters */
lhausamm's avatar
lhausamm committed
191
  parser_get_param_string(parameter_file, "GrackleCooling:GrackleCloudyTable",
192
193
194
195
196
197
198
199
                          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
200
      parameter_file, "GrackleCooling:GrackleHSShieldingDensityThreshold");
201
202
203
204
205
  
#ifdef SWIFT_DEBUG_CHECKS
  /* enable verbose for grackle */
  grackle_verbose = 1;
#endif
lhausamm's avatar
lhausamm committed
206

207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
  /* 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)  */
  cooling->comoving_coordinates = 0;

  /* 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;
 
  /* Create a chemistry object for parameters and rate data. */
  if (set_default_chemistry_parameters() == 0) {
    error("Error in set_default_chemistry_parameters.");
  }
lhausamm's avatar
lhausamm committed
228

229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
  // Set parameter values for chemistry.
  grackle_data.use_grackle = 1;
  grackle_data.with_radiative_cooling = 1;
  /* molecular network with H, He, D
   From Cloudy table */
  grackle_data.primordial_chemistry = 0;
  grackle_data.metal_cooling = 1;  // metal cooling on
  grackle_data.UVbackground = cooling->uv_background;
  grackle_data.grackle_data_file = cooling->cloudy_table;

  /* Initialize the chemistry object.
     a_value is not the true initial a
     This should get set before any computation */
  double a_value = 1.;

  if (initialize_chemistry_data(&cooling->units, a_value) == 0) {
    error("Error in initialize_chemistry_data.");
  }
lhausamm's avatar
lhausamm committed
247
248

#ifdef SWIFT_DEBUG_CHECKS
249
250
  if (GRACKLE_NPART != 1)
    error("Grackle with multiple particles not implemented");
lhausamm's avatar
lhausamm committed
251
252
253
254
255
256
257
258
  float threshold = cooling->GrackleHSShieldingDensityThreshold;

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

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

262

263
  //grackle_print_data();
lhausamm's avatar
lhausamm committed
264
265
266
  message("");
  message("***************************************");
#endif
267
268
269
270
271
272
273
274
275
276
277
}

/**
 * @brief Prints the properties of the cooling model to stdout.
 *
 * @param cooling The properties of the cooling function.
 */
static INLINE void cooling_print_backend(
    const struct cooling_function_data* cooling) {

  message("Cooling function is 'Grackle'.");
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
  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     = %g", cooling->units.comoving_coordinates)
  message("\tLength  = %g", cooling->units.length_units);
  message("\tDensity = %g", cooling->units.density_units);
  message("\tTime         = %g", cooling->units.time_units);
  message("\tScale Factor = %g", cooling->units.a_units);  
}


/**
 * @brief print data in grackle struct
 *
 * Should only be used for debugging
 */
void grackle_print_data() {
  message("Grackle Data:");
  message("\t Data file: %s", grackle_data.grackle_data_file);
  message("\t With grackle: %i", grackle_data.use_grackle);
  message("\t With radiative cooling: %i", grackle_data.with_radiative_cooling);
  message("\t With UV background: %i", grackle_data.UVbackground);
  message("\t With primordial chemistry: %i",
          grackle_data.primordial_chemistry);
  message("\t Number temperature bins: %i",
          grackle_data.NumberOfTemperatureBins);
  message("\t T = (%g, ..., %g)", grackle_data.TemperatureStart,
          grackle_data.TemperatureEnd);

  message("Primordial Cloudy");
  cloudy_print_data(grackle_data.cloudy_primordial, 1);
  if (grackle_data.metal_cooling) {
    message("Metal Cooling");
    cloudy_print_data(grackle_data.cloudy_metal, 0);
  }

  message("\t Gamma: %g", grackle_data.Gamma);

  /* UVB */
  if (grackle_data.UVbackground && grackle_data.primordial_chemistry != 0) {
    struct UVBtable uvb = grackle_data.UVbackground_table;
    long long N = uvb.Nz;
    message("\t UV Background");
    message("\t\t Redshift from %g to %g with %lli steps", uvb.zmin, uvb.zmax,
            N);
    message("\t\t z = (%g, ..., %g)", uvb.z[0], uvb.z[N - 1]);
  }
}

/**
 * @brief print data in cloudy struct
 *
 * Should only be used for debugging
 */
void cloudy_print_data(const cloudy_data c, const int print_mmw) {
  long long N = c.data_size;
  message("\t Data size: %lli", N);
  message("\t Grid rank: %lli", c.grid_rank);

  char msg[200] = "\t Dimension: (";
  for (long long i = 0; i < c.grid_rank; i++) {
    char tmp[200] = "%lli%s";
    if (i == c.grid_rank - 1)
      sprintf(tmp, tmp, c.grid_dimension[i], ")");
    else
      sprintf(tmp, tmp, c.grid_dimension[i], ", ");

    strcat(msg, tmp);
  }
  message("%s", msg);

  if (c.heating_data)
    message("\t Heating: (%g, ..., %g)", c.heating_data[0],
            c.heating_data[N - 1]);
  if (c.cooling_data)
    message("\t Cooling: (%g, ..., %g)", c.cooling_data[0],
            c.cooling_data[N - 1]);
  if (c.mmw_data && print_mmw)
    message("\t Mean molecular weigth: (%g, ..., %g)", c.mmw_data[0],
            c.mmw_data[N - 1]);
362
363
364
}

#endif /* SWIFT_COOLING_GRACKLE_H */