cooling.h 21.6 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

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

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

46
47
#ifdef HAVE_HDF5

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

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

59
/**
60
61
 * @brief Sets the cooling properties of the (x-)particles to a valid start
 * state.
62
 *
63
64
 * @param p Pointer to the particle data.
 * @param xp Pointer to the extended particle data.
65
 * @param cooling The properties of the cooling function.
66
 */
67
68
69
__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) {
70
71

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

73
74
75
  /* metal cooling = 1 */
  xp->cooling_data.metal_frac = cooling->chemistry.SolarMetalFractionByMass;

lhausamm's avatar
lhausamm committed
76
#if COOLING_GRACKLE_MODE >= 1
lhausamm's avatar
lhausamm committed
77
78
  gr_float zero = 1.e-20;

lhausamm's avatar
lhausamm committed
79
  /* primordial chemistry >= 1 */
80
  xp->cooling_data.HI_frac = grackle_data->HydrogenFractionByMass;
lhausamm's avatar
lhausamm committed
81
  xp->cooling_data.HII_frac = zero;
82
  xp->cooling_data.HeI_frac = 1. - grackle_data->HydrogenFractionByMass;
lhausamm's avatar
lhausamm committed
83
84
85
  xp->cooling_data.HeII_frac = zero;
  xp->cooling_data.HeIII_frac = zero;
  xp->cooling_data.e_frac = zero;
lhausamm's avatar
lhausamm committed
86
87
88

#if COOLING_GRACKLE_MODE >= 2
  /* primordial chemistry >= 2 */
lhausamm's avatar
lhausamm committed
89
90
91
  xp->cooling_data.HM_frac = zero;
  xp->cooling_data.H2I_frac = zero;
  xp->cooling_data.H2II_frac = zero;
lhausamm's avatar
lhausamm committed
92
93
94

#if COOLING_GRACKLE_MODE >= 3
  /* primordial chemistry >= 3 */
lhausamm's avatar
lhausamm committed
95
96
  xp->cooling_data.DI_frac = grackle_data->HydrogenFractionByMass *
                             grackle_data->DeuteriumToHydrogenRatio;
lhausamm's avatar
lhausamm committed
97
98
  xp->cooling_data.DII_frac = zero;
  xp->cooling_data.HDI_frac = zero;
lhausamm's avatar
lhausamm committed
99
#endif  // MODE >= 3
lhausamm's avatar
lhausamm committed
100

lhausamm's avatar
lhausamm committed
101
#endif  // MODE >= 2
lhausamm's avatar
lhausamm committed
102

lhausamm's avatar
lhausamm committed
103
#endif  // MODE >= 1
104
105
106
107
108
109
110
111
}

/**
 * @brief update particle with densities.
 *
 * @param xp The extended particle data
 */
__attribute__((always_inline)) INLINE static void cooling_compute_density(
lhausamm's avatar
lhausamm committed
112
    struct xpart* restrict xp, const gr_float rho) {
113
114
115
116
117
118
119
120
121
122

#if COOLING_GRACKLE_MODE >= 1
  /* primordial chemistry >= 1 */
  xp->cooling_data.HI_frac *= rho;
  xp->cooling_data.HII_frac *= rho;
  xp->cooling_data.HeI_frac *= rho;
  xp->cooling_data.HeII_frac *= rho;
  xp->cooling_data.HeIII_frac *= rho;
  xp->cooling_data.e_frac *= rho;

lhausamm's avatar
lhausamm committed
123
#if COOLING_GRACKLE_MODE >= 2
124
125
126
127
128
129
130
131
132
133
  /* primordial chemistry >= 2 */
  xp->cooling_data.HM_frac *= rho;
  xp->cooling_data.H2I_frac *= rho;
  xp->cooling_data.H2II_frac *= rho;

#if COOLING_GRACKLE_MODE >= 3
  /* primordial chemistry >= 3 */
  xp->cooling_data.DI_frac *= rho;
  xp->cooling_data.DII_frac *= rho;
  xp->cooling_data.HDI_frac *= rho;
lhausamm's avatar
lhausamm committed
134
#endif  // MODE >= 3
135

lhausamm's avatar
lhausamm committed
136
#endif  // MODE >= 2
137

lhausamm's avatar
lhausamm committed
138
#endif  // MODE >= 1
139
140
141
142
143
144
145
146
147
148

  xp->cooling_data.metal_frac *= rho;
}

/**
 * @brief update particle with fraction.
 *
 * @param xp The extended particle data
 */
__attribute__((always_inline)) INLINE static void cooling_compute_fraction(
lhausamm's avatar
lhausamm committed
149
    struct xpart* restrict xp, const gr_float rho) {
150
151
152
153
154
155
156
157
158
159

#if COOLING_GRACKLE_MODE >= 1
  /* primordial chemistry >= 1 */
  xp->cooling_data.HI_frac /= rho;
  xp->cooling_data.HII_frac /= rho;
  xp->cooling_data.HeI_frac /= rho;
  xp->cooling_data.HeII_frac /= rho;
  xp->cooling_data.HeIII_frac /= rho;
  xp->cooling_data.e_frac /= rho;

lhausamm's avatar
lhausamm committed
160
#if COOLING_GRACKLE_MODE >= 2
161
162
163
164
165
166
167
168
169
170
  /* primordial chemistry >= 2 */
  xp->cooling_data.HM_frac /= rho;
  xp->cooling_data.H2I_frac /= rho;
  xp->cooling_data.H2II_frac /= rho;

#if COOLING_GRACKLE_MODE >= 3
  /* primordial chemistry >= 3 */
  xp->cooling_data.DI_frac /= rho;
  xp->cooling_data.DII_frac /= rho;
  xp->cooling_data.HDI_frac /= rho;
lhausamm's avatar
lhausamm committed
171
#endif  // MODE >= 3
172

lhausamm's avatar
lhausamm committed
173
#endif  // MODE >= 2
174

lhausamm's avatar
lhausamm committed
175
#endif  // MODE >= 1
176
177

  xp->cooling_data.metal_frac /= rho;
178
179
180
181
182
183
184
185
186
187
188
189
190
}

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

191
192
193
194
195
/**
 * @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
196
__attribute__((always_inline)) INLINE static void cooling_print_backend(
197
198
199
    const struct cooling_function_data* cooling) {

  message("Cooling function is 'Grackle'.");
200
  message("Using Grackle           = %i", cooling->chemistry.use_grackle);
Matthieu Schaller's avatar
Matthieu Schaller committed
201
202
203
204
  message("Chemical network        = %i",
          cooling->chemistry.primordial_chemistry);
  message("Radiative cooling       = %i",
          cooling->chemistry.with_radiative_cooling);
205
  message("Metal cooling           = %i", cooling->chemistry.metal_cooling);
Matthieu Schaller's avatar
Matthieu Schaller committed
206
207

  message("CloudyTable             = %s", cooling->cloudy_table);
208
209
  message("UVbackground            = %d", cooling->uv_background);
  message("Redshift                = %g", cooling->redshift);
lhausamm's avatar
lhausamm committed
210
211
  message("Solar Metal Fraction    = %g",
          cooling->chemistry.SolarMetalFractionByMass);
212
213
  message("Units:");
  message("\tComoving     = %i", cooling->units.comoving_coordinates);
lhausamm's avatar
lhausamm committed
214
215
  message("\tLength       = %g", cooling->units.length_units);
  message("\tDensity      = %g", cooling->units.density_units);
216
  message("\tTime         = %g", cooling->units.time_units);
lhausamm's avatar
lhausamm committed
217
  message("\tScale Factor = %g", cooling->units.a_units);
218
#ifdef SWIFT_DEBUG_CHECKS
lhausamm's avatar
lhausamm committed
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
/*
const chemistry_data *tmp = &cooling->chemistry;
message("Debug:");
message("UVBackground                       = %i", tmp->UVbackground);
message("Grackle data file                  = %s", tmp->grackle_data_file);
message("CMB temperature floor              = %i", tmp->cmb_temperature_floor);
message("Gamma                              = %g", tmp->Gamma);
message("H2 on dust                         = %i", tmp->h2_on_dust);
message("Photoelectric heating              = %i", tmp->photoelectric_heating);
message("Photoelectric heating rate         = %g",
tmp->photoelectric_heating_rate);
message("Use volumetric heating rate        = %i",
tmp->use_volumetric_heating_rate);
message("Use specific heating rate          = %i",
tmp->use_specific_heating_rate);
message("Three body                         = %i", tmp->three_body_rate);
message("Cie cooling                        = %i", tmp->cie_cooling);
message("h2 optical depth approx            = %i",
tmp->h2_optical_depth_approximation);
message("ih2co                              = %i", tmp->ih2co);
message("ipiht                              = %i", tmp->ipiht);

message("Hydrogen Fraction                  = %g", tmp->HydrogenFractionByMass);
message("Deuterium/Hydrogen ratio           = %g",
tmp->DeuteriumToHydrogenRatio);
message("Solar metal fraction               = %g",
tmp->SolarMetalFractionByMass);

message("Number T bins                      = %i",
tmp->NumberOfTemperatureBins);
message("Case B recombination               = %i", tmp->CaseBRecombination);

message("T start                            = %g", tmp->TemperatureStart);
message("T end                              = %g", tmp->TemperatureEnd);

message("Number dust T bins                 = %i",
tmp->NumberOfDustTemperatureBins);
message("Dust T start                       = %g", tmp->DustTemperatureStart);
message("Dust T end                         = %g", tmp->DustTemperatureEnd);

message("Compton xray heating               = %i", tmp->Compton_xray_heating);
message("LW background sawtooth suppression = %i",
tmp->LWbackground_sawtooth_suppression);
message("LW background intensity            = %g", tmp->LWbackground_intensity);
message("UV redshift on                     = %g",
tmp->UVbackground_redshift_on);
message("UV redshift off                    = %g",
tmp->UVbackground_redshift_off);
message("UV redshift fullon                 = %g",
tmp->UVbackground_redshift_fullon);
message("UV redshift drop                   = %g",
tmp->UVbackground_redshift_drop);

message("Cloudy electron fraction           = %g",
tmp->cloudy_electron_fraction_factor);

message("Use radiative transfer             = %i", tmp->use_radiative_transfer);
message("RT coupled rate solver             = %i",
tmp->radiative_transfer_coupled_rate_solver);
message("RT intermediate step               = %i",
tmp->radiative_transfer_intermediate_step);
message("RT H only                          = %i",
tmp->radiative_transfer_hydrogen_only);

message("Self shielding method              = %i", tmp->self_shielding_method);
*/
285
#endif
286
287
}

lhausamm's avatar
lhausamm committed
288
289
290
291
292
/**
 * @brief allocate required field
 *
 * @param data the #grackle_field_data
 */
lhausamm's avatar
lhausamm committed
293
294
__attribute__((always_inline)) INLINE static void cooling_malloc_data(
    grackle_field_data* data) {
lhausamm's avatar
lhausamm committed
295
296
297
298
299
300
301
302
303

#if COOLING_GRACKLE_MODE >= 1
  /* primordial chemistry >= 1 */
  data->HI_density = malloc(sizeof(gr_float));
  data->HII_density = malloc(sizeof(gr_float));
  data->HeI_density = malloc(sizeof(gr_float));
  data->HeII_density = malloc(sizeof(gr_float));
  data->HeIII_density = malloc(sizeof(gr_float));
  data->e_density = malloc(sizeof(gr_float));
lhausamm's avatar
lhausamm committed
304
#endif  // MODE >= 1
lhausamm's avatar
lhausamm committed
305
306
307
308
309
310

#if COOLING_GRACKLE_MODE >= 2
  /* primordial chemistry >= 2 */
  data->HM_density = malloc(sizeof(gr_float));
  data->H2I_density = malloc(sizeof(gr_float));
  data->H2II_density = malloc(sizeof(gr_float));
lhausamm's avatar
lhausamm committed
311
#endif  // MODE 2
lhausamm's avatar
lhausamm committed
312
313
314
315
316
317

#if COOLING_GRACKLE_MODE >= 3
  /* primordial chemistry >= 3 */
  data->DI_density = malloc(sizeof(gr_float));
  data->DII_density = malloc(sizeof(gr_float));
  data->HDI_density = malloc(sizeof(gr_float));
lhausamm's avatar
lhausamm committed
318
#endif  // MODE >= 3
lhausamm's avatar
lhausamm committed
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335

  /* metal cooling = 1 */
  data->metal_density = malloc(sizeof(gr_float));

  /* /\* volumetric heating rate *\/ */
  /* data->volumetric_heating_rate = NULL; */

  /* /\* specific heating rate *\/ */
  /* data->specific_heating_rate = NULL; */
}

/**
 * @brief free the allocated memory
 *
 * @param data the #grackle_field_data
 */

lhausamm's avatar
lhausamm committed
336
337
__attribute__((always_inline)) INLINE static void cooling_free_data(
    grackle_field_data* data) {
lhausamm's avatar
lhausamm committed
338
339
340
341
342
343
344
345
346

#if COOLING_GRACKLE_MODE >= 1
  /* primordial chemistry >= 1 */
  free(data->HI_density);
  free(data->HII_density);
  free(data->HeI_density);
  free(data->HeII_density);
  free(data->HeIII_density);
  free(data->e_density);
lhausamm's avatar
lhausamm committed
347
#endif  // MODE >= 1
lhausamm's avatar
lhausamm committed
348
349
350
351
352
353

#if COOLING_GRACKLE_MODE >= 2
  /* primordial chemistry >= 2 */
  free(data->HM_density);
  free(data->H2I_density);
  free(data->H2II_density);
lhausamm's avatar
lhausamm committed
354
#endif  // MODE 2
lhausamm's avatar
lhausamm committed
355
356
357
358
359
360

#if COOLING_GRACKLE_MODE >= 3
  /* primordial chemistry >= 3 */
  free(data->DI_density);
  free(data->DII_density);
  free(data->HDI_density);
lhausamm's avatar
lhausamm committed
361
#endif  // MODE >= 3
lhausamm's avatar
lhausamm committed
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380

  /* metal cooling = 1 */
  free(data->metal_density);

  /* /\* volumetric heating rate *\/ */
  /* data->volumetric_heating_rate = NULL; */

  /* /\* specific heating rate *\/ */
  /* data->specific_heating_rate = NULL; */
}

/**
 * @brief copy xp to data
 *
 * requires the particle to have been transformed into density
 *
 * @param data the #grackle_field_data
 * @param xp the #xpart
 */
lhausamm's avatar
lhausamm committed
381
382
__attribute__((always_inline)) INLINE static void cooling_copy_to_data(
    grackle_field_data* data, const struct xpart* xp) {
lhausamm's avatar
lhausamm committed
383
384
385
386
387
388
389
390
391

#if COOLING_GRACKLE_MODE >= 1
  /* primordial chemistry >= 1 */
  data->HI_density[0] = xp->cooling_data.HI_frac;
  data->HII_density[0] = xp->cooling_data.HII_frac;
  data->HeI_density[0] = xp->cooling_data.HeI_frac;
  data->HeII_density[0] = xp->cooling_data.HeII_frac;
  data->HeIII_density[0] = xp->cooling_data.HeIII_frac;
  data->e_density[0] = xp->cooling_data.e_frac;
lhausamm's avatar
lhausamm committed
392
#endif  // MODE >= 1
lhausamm's avatar
lhausamm committed
393
394
395
396
397
398

#if COOLING_GRACKLE_MODE >= 2
  /* primordial chemistry >= 2 */
  data->HM_density[0] = xp->cooling_data.HM_frac;
  data->H2I_density[0] = xp->cooling_data.H2I_frac;
  data->H2II_density[0] = xp->cooling_data.H2II_frac;
lhausamm's avatar
lhausamm committed
399
#endif  // MODE 2
lhausamm's avatar
lhausamm committed
400
401
402
403
404
405

#if COOLING_GRACKLE_MODE >= 3
  /* primordial chemistry >= 3 */
  data->DI_density[0] = xp->cooling_data.DI_frac;
  data->DII_density[0] = xp->cooling_data.DII_frac;
  data->HDI_density[0] = xp->cooling_data.HDI_frac;
lhausamm's avatar
lhausamm committed
406
#endif  // MODE >= 3
lhausamm's avatar
lhausamm committed
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423

  /* metal cooling = 1 */
  data->metal_density[0] = xp->cooling_data.metal_frac;

  /* volumetric heating rate */
  data->volumetric_heating_rate = NULL;

  /* specific heating rate */
  data->specific_heating_rate = NULL;
}

/**
 * @brief copy data to xp
 *
 * @param data the #grackle_field_data
 * @param xp the #xpart
 */
lhausamm's avatar
lhausamm committed
424
425
__attribute__((always_inline)) INLINE static void cooling_copy_to_particle(
    const grackle_field_data* data, struct xpart* xp) {
lhausamm's avatar
lhausamm committed
426
427
428
429
430
431
432
433
434

#if COOLING_GRACKLE_MODE >= 1
  /* primordial chemistry >= 1 */
  xp->cooling_data.HI_frac = data->HI_density[0];
  xp->cooling_data.HII_frac = data->HII_density[0];
  xp->cooling_data.HeI_frac = data->HeI_density[0];
  xp->cooling_data.HeII_frac = data->HeII_density[0];
  xp->cooling_data.HeIII_frac = data->HeIII_density[0];
  xp->cooling_data.e_frac = data->e_density[0];
lhausamm's avatar
lhausamm committed
435
#endif  // MODE >= 1
lhausamm's avatar
lhausamm committed
436
437
438
439
440
441

#if COOLING_GRACKLE_MODE >= 2
  /* primordial chemistry >= 2 */
  xp->cooling_data.HM_frac = data->HM_density[0];
  xp->cooling_data.H2I_frac = data->H2I_density[0];
  xp->cooling_data.H2II_frac = data->H2II_density[0];
lhausamm's avatar
lhausamm committed
442
#endif  // MODE 2
lhausamm's avatar
lhausamm committed
443
444
445
446
447
448

#if COOLING_GRACKLE_MODE >= 3
  /* primordial chemistry >= 3 */
  xp->cooling_data.DI_frac = data->DI_density[0];
  xp->cooling_data.DII_frac = data->DII_density[0];
  xp->cooling_data.HDI_frac = data->HDI_density[0];
lhausamm's avatar
lhausamm committed
449
#endif  // MODE >= 3
lhausamm's avatar
lhausamm committed
450
451
452
453
454
455
456
457
458
459
460

  /* metal cooling = 1 */
  xp->cooling_data.metal_frac = data->metal_density[0];

  /* /\* volumetric heating rate *\/ */
  /* data->volumetric_heating_rate = NULL; */

  /* /\* specific heating rate *\/ */
  /* data->specific_heating_rate = NULL; */
}

461
/**
462
 * @brief Compute the cooling rate and update the particle chemistry data
463
464
465
466
467
 *
 * @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.
468
 * @param xp Pointer to the particle extra data
469
 * @param dt The time-step of this particle.
470
471
 *
 * @return du / dt
472
 */
473
474
475
__attribute__((always_inline)) INLINE static double cooling_rate(
    const struct phys_const* restrict phys_const,
    const struct unit_system* restrict us,
476
    const struct cosmology* restrict cosmo,
477
    const struct cooling_function_data* restrict cooling,
lhausamm's avatar
lhausamm committed
478
    struct part* restrict p, struct xpart* restrict xp, double dt) {
479

lhausamm's avatar
lhausamm committed
480
  /* set current time */
lhausamm's avatar
lhausamm committed
481
  code_units units = cooling->units;
lhausamm's avatar
lhausamm committed
482
483
484
  if (cooling->redshift == -1)
    error("TODO time dependant redshift");
  else
lhausamm's avatar
lhausamm committed
485
    units.a_value = 1. / (1. + cooling->redshift);
486

lhausamm's avatar
lhausamm committed
487
488
  /* initialize data */
  grackle_field_data data;
lhausamm's avatar
lhausamm committed
489

lhausamm's avatar
lhausamm committed
490
491
492
493
494
  /* 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
495

lhausamm's avatar
lhausamm committed
496
497
498
499
  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
500

lhausamm's avatar
lhausamm committed
501
  /* general particle data */
502
503
  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
504
  gr_float energy = energy_before;
505

lhausamm's avatar
lhausamm committed
506
  /* initialize density */
lhausamm's avatar
lhausamm committed
507
  data.density = &density;
lhausamm's avatar
lhausamm committed
508
509

  /* initialize energy */
lhausamm's avatar
lhausamm committed
510
  data.internal_energy = &energy;
lhausamm's avatar
lhausamm committed
511
512
513
514
515

  /* grackle 3.0 doc: "Currently not used" */
  data.x_velocity = NULL;
  data.y_velocity = NULL;
  data.z_velocity = NULL;
lhausamm's avatar
lhausamm committed
516

517
  /* transform gas fraction to densities */
lhausamm's avatar
lhausamm committed
518
  cooling_compute_density(xp, *data.density);
519

lhausamm's avatar
lhausamm committed
520
521
  /* allocate grackle data */
  cooling_malloc_data(&data);
522

lhausamm's avatar
lhausamm committed
523
524
  /* copy data from particle to grackle data */
  cooling_copy_to_data(&data, xp);
lhausamm's avatar
lhausamm committed
525
526

  /* solve chemistry with table */
lhausamm's avatar
lhausamm committed
527
  if (solve_chemistry(&units, &data, dt) == 0) {
lhausamm's avatar
lhausamm committed
528
    error("Error in solve_chemistry.");
529
  }
lhausamm's avatar
lhausamm committed
530

lhausamm's avatar
lhausamm committed
531
532
533
  /* copy from grackle data to particle */
  cooling_copy_to_particle(&data, xp);

534
  /* transform densities to gas fraction */
lhausamm's avatar
lhausamm committed
535
  cooling_compute_fraction(xp, *data.density);
536

lhausamm's avatar
lhausamm committed
537
538
539
  /* free allocated memory */
  cooling_free_data(&data);

lhausamm's avatar
lhausamm committed
540
  /* compute rate */
lhausamm's avatar
lhausamm committed
541
  return (energy - energy_before) / dt;
542
543
544
545
546
547
548
}

/**
 * @brief Apply the cooling function to a particle.
 *
 * @param phys_const The physical constants in internal units.
 * @param us The internal system of units.
549
 * @param cosmo The current cosmological model.
550
551
552
553
554
555
 * @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
556
    const struct unit_system* restrict us,
557
    const struct cosmology* restrict cosmo,
558
    const struct cooling_function_data* restrict cooling,
lhausamm's avatar
lhausamm committed
559
    struct part* restrict p, struct xpart* restrict xp, double dt) {
560

Matthieu Schaller's avatar
Matthieu Schaller committed
561
  if (dt == 0.) return;
562
563
564

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

566
  /* compute cooling rate */
567
  const float du_dt = cooling_rate(phys_const, us, cosmo, cooling, p, dt);
568
569

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

572
  /* Update the internal energy */
573
  hydro_set_internal_energy_dt(p, hydro_du_dt + du_dt);
574
575
576
577
578
579
580
581
582
}

/**
 * @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.
583
 * @param cosmo The current cosmological model.
584
585
586
587
588
589
 * @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,
590
    const struct cosmology* restrict cosmo,
lhausamm's avatar
lhausamm committed
591
    const struct unit_system* restrict us, const struct part* restrict p) {
592
593
594
595
596
597
598
599
600
601
602
603

  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
604
__attribute__((always_inline)) INLINE static void cooling_init_backend(
lhausamm's avatar
lhausamm committed
605
    const struct swift_params* parameter_file, const struct unit_system* us,
606
607
608
    const struct phys_const* phys_const,
    struct cooling_function_data* cooling) {

lhausamm's avatar
lhausamm committed
609
610
611
  if (GRACKLE_NPART != 1)
    error("Grackle with multiple particles not implemented");

Matthieu Schaller's avatar
Matthieu Schaller committed
612
  /* read parameters */
lhausamm's avatar
lhausamm committed
613
  parser_get_param_string(parameter_file, "GrackleCooling:GrackleCloudyTable",
614
615
                          cooling->cloudy_table);
  cooling->uv_background =
Matthieu Schaller's avatar
Matthieu Schaller committed
616
      parser_get_param_int(parameter_file, "GrackleCooling:UVbackground");
617
618

  cooling->redshift =
Matthieu Schaller's avatar
Matthieu Schaller committed
619
      parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift");
620
621
622
623
624
625
626
627
628
629

#ifdef SWIFT_DEBUG_CHECKS
  /* enable verbose for grackle */
  grackle_verbose = 1;
#endif
  /* 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)
630
  cooling->units.a_value = 1.0;
631
632
633

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

  /* then units */
Matthieu Schaller's avatar
Matthieu Schaller committed
637
638
  cooling->units.density_units =
      us->UnitMass_in_cgs / pow(us->UnitLength_in_cgs, 3);
639
640
  cooling->units.length_units = us->UnitLength_in_cgs;
  cooling->units.time_units = us->UnitTime_in_cgs;
Matthieu Schaller's avatar
Matthieu Schaller committed
641
642
643
  cooling->units.velocity_units = cooling->units.a_units *
                                  cooling->units.length_units /
                                  cooling->units.time_units;
lhausamm's avatar
lhausamm committed
644

Matthieu Schaller's avatar
Matthieu Schaller committed
645
  chemistry_data* chemistry = &cooling->chemistry;
lhausamm's avatar
lhausamm committed
646

647
  /* Create a chemistry object for parameters and rate data. */
lhausamm's avatar
lhausamm committed
648
  if (set_default_chemistry_parameters(chemistry) == 0) {
649
650
    error("Error in set_default_chemistry_parameters.");
  }
lhausamm's avatar
lhausamm committed
651

652
  // Set parameter values for chemistry.
lhausamm's avatar
lhausamm committed
653
654
  chemistry->use_grackle = 1;
  chemistry->with_radiative_cooling = 1;
lhausamm's avatar
lhausamm committed
655

656
657
  /* molecular network with H, He, D
   From Cloudy table */
658
  chemistry->primordial_chemistry = COOLING_GRACKLE_MODE;
659
  chemistry->metal_cooling = 1;
lhausamm's avatar
lhausamm committed
660
661
  chemistry->UVbackground = cooling->uv_background;
  chemistry->grackle_data_file = cooling->cloudy_table;
lhausamm's avatar
lhausamm committed
662

663
664
665
  chemistry->use_radiative_transfer = 0;
  chemistry->use_volumetric_heating_rate = 0;
  chemistry->use_specific_heating_rate = 0;
lhausamm's avatar
lhausamm committed
666

lhausamm's avatar
lhausamm committed
667
668
  /* Initialize the chemistry object. */
  if (initialize_chemistry_data(&cooling->units) == 0) {
669
670
    error("Error in initialize_chemistry_data.");
  }
lhausamm's avatar
lhausamm committed
671
672
673
674
675

#ifdef SWIFT_DEBUG_CHECKS
  message("***************************************");
  message("initializing grackle cooling function");
  message("");
676
  cooling_print_backend(cooling);
lhausamm's avatar
lhausamm committed
677
678
679
  message("");
  message("***************************************");
#endif
680
681
682
}

#endif /* SWIFT_COOLING_GRACKLE_H */