cooling.h 24.8 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
#include "chemistry.h"
lhausamm's avatar
lhausamm committed
35
#include "cooling_io.h"
36
37
#include "error.h"
#include "hydro.h"
lhausamm's avatar
lhausamm committed
38
#include "io_properties.h"
39
40
41
42
43
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
#include "units.h"

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

48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
/* prototypes */
static gr_float cooling_time(
    const struct cooling_function_data* restrict cooling,
    const struct part* restrict p, struct xpart* restrict xp);

static double cooling_rate(
    const struct phys_const* restrict phys_const,
    const struct unit_system* restrict us,
    const struct cooling_function_data* restrict cooling,
    const struct part* restrict p, struct xpart* restrict xp, double dt);

/**
 * @brief Print the chemical network
 *
 * @param xp The #xpart to print
 */
__attribute__((always_inline)) INLINE static void cooling_print_fractions(
65
    const struct xpart* restrict xp) {
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

  const struct cooling_xpart_data tmp = xp->cooling_data;
#if COOLING_GRACKLE_MODE > 0
  message("HI %g, HII %g, HeI %g, HeII %g, HeIII %g, e %g",
	  tmp.HI_frac, tmp.HII_frac, tmp.HeI_frac,
	  tmp.HeII_frac, tmp.HeIII_frac, tmp.e_frac);
#endif

#if COOLING_GRACKLE_MODE > 1
  message("HM %g, H2I %g, H2II %g",
	  tmp.HM_frac, tmp.H2I_frac, tmp.H2II_frac);
#endif

#if COOLING_GRACKLE_MODE > 2
  message("DI %g, DII %g, HDI %g",
	  tmp.DI_frac, tmp.DII_frac, tmp.HDI_frac);
#endif
  message("Metal: %g", tmp.metal_frac);
}

/**
lhausamm's avatar
lhausamm committed
87
 * @brief Check if the difference of a given field is lower than limit
88
 *
lhausamm's avatar
lhausamm committed
89
90
91
92
 * @param xp First #xpart
 * @param old Second #xpart
 * @param field The field to check
 * @param limit Difference limit
93
 *
lhausamm's avatar
lhausamm committed
94
 * @return 0 if diff > limit
95
 */
lhausamm's avatar
lhausamm committed
96
97
98
99
100
101
102
#define cooling_check_field(xp, old, field, limit)			\
  ({									\
    float tmp = xp->cooling_data.field - old->cooling_data.field;	\
    tmp = fabs(tmp) / xp->cooling_data.field;				\
    if (tmp > limit)							\
      return 0;								\
})
103

lhausamm's avatar
lhausamm committed
104
105
106
107
108
109
110
111
112
113
114
/**
 * @brief Check if difference between two particles is lower than a given value
 *
 * @param xp One of the #xpart
 * @param old The other #xpart
 * @param limit The difference limit
 */
__attribute__((always_inline)) INLINE static int cooling_converged(
    const struct xpart* restrict xp,
    const struct xpart* restrict old,
    const float limit) {
115
116

#if COOLING_GRACKLE_MODE > 0
lhausamm's avatar
lhausamm committed
117
118
119
120
121
122
123
  cooling_check_field(xp, old, HI_frac, limit);
  cooling_check_field(xp, old, HII_frac, limit);
  cooling_check_field(xp, old, HeI_frac, limit);
  cooling_check_field(xp, old, HeII_frac, limit);
  cooling_check_field(xp, old, HeIII_frac, limit);
  cooling_check_field(xp, old, e_frac, limit);
#endif
124
#if COOLING_GRACKLE_MODE > 1
lhausamm's avatar
lhausamm committed
125
126
127
128
  cooling_check_field(xp, old, HM_frac, limit);
  cooling_check_field(xp, old, H2I_frac, limit);
  cooling_check_field(xp, old, H2II_frac, limit);
#endif
129
130

#if COOLING_GRACKLE_MODE > 2
lhausamm's avatar
lhausamm committed
131
132
133
134
  cooling_check_field(xp, old, DI_frac, limit);
  cooling_check_field(xp, old, DII_frac, limit);
  cooling_check_field(xp, old, HDI_frac, limit);
#endif
135

136
  return 1;
137
138
139
}

/**
lhausamm's avatar
lhausamm committed
140
 * @brief Compute equilibrium fraction
141
 *
lhausamm's avatar
lhausamm committed
142
143
144
 * @param p Pointer to the particle data.
 * @param xp Pointer to the extended particle data.
 * @param cooling The properties of the cooling function.
145
 */
lhausamm's avatar
lhausamm committed
146
147
148
149
__attribute__((always_inline)) INLINE static void cooling_compute_equilibrium(
    const struct part* restrict p,
    struct xpart* restrict xp,
    const struct cooling_function_data* cooling) {
150

lhausamm's avatar
lhausamm committed
151
152
153
154
155
156
  /* get temporary data */
  struct part p_tmp = *p;
  struct cooling_function_data cooling_tmp = *cooling;
  cooling_tmp.chemistry.with_radiative_cooling = 0;
  /* need density for computation, therefore quick estimate */
  p_tmp.rho = 0.2387 * p_tmp.mass / pow(p_tmp.h, 3);
157

lhausamm's avatar
lhausamm committed
158
159
160
161
162
  /* compute time step */
  const double alpha = 0.01;
  double dt = fabs(cooling_time(&cooling_tmp, &p_tmp, xp));
  cooling_rate(NULL, NULL, &cooling_tmp, &p_tmp, xp, dt);
  dt = alpha * fabs(cooling_time(&cooling_tmp, &p_tmp, xp));
163

lhausamm's avatar
lhausamm committed
164
  /* init simple variables */
165
  int step = 0;
lhausamm's avatar
lhausamm committed
166
167
168
  const int max_step = cooling_tmp.max_step;
  const float conv_limit = cooling_tmp.convergence_limit;
  struct xpart old;
169
170

  do {
lhausamm's avatar
lhausamm committed
171
    /* update variables */
172
    step += 1;
lhausamm's avatar
lhausamm committed
173
    old = *xp;
174

lhausamm's avatar
lhausamm committed
175
176
177
    /* update chemistry */
    cooling_rate(NULL, NULL, &cooling_tmp, &p_tmp, xp, dt);
  } while (step < max_step && !cooling_converged(xp, &old, conv_limit));
178

lhausamm's avatar
lhausamm committed
179
180
181
182
  if (step == max_step)
    error("A particle element fraction failed to converge."
	  "You can change 'GrackleCooling:MaxSteps' or "
	  "'GrackleCooling:ConvergenceLimit' to avoid this problem");
183
184
}

185
/**
186
187
 * @brief Sets the cooling properties of the (x-)particles to a valid start
 * state.
188
 *
189
190
 * @param p Pointer to the particle data.
 * @param xp Pointer to the extended particle data.
191
 * @param cooling The properties of the cooling function.
192
 */
193
194
195
__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) {
196
197

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

lhausamm's avatar
lhausamm committed
199
#if COOLING_GRACKLE_MODE >= 1
lhausamm's avatar
lhausamm committed
200
201
  gr_float zero = 1.e-20;

lhausamm's avatar
lhausamm committed
202
  /* primordial chemistry >= 1 */
203
204
  xp->cooling_data.HI_frac = zero;
  xp->cooling_data.HII_frac = grackle_data->HydrogenFractionByMass;
205
  xp->cooling_data.HeI_frac = 1. - grackle_data->HydrogenFractionByMass;
lhausamm's avatar
lhausamm committed
206
207
  xp->cooling_data.HeII_frac = zero;
  xp->cooling_data.HeIII_frac = zero;
208
209
210
211
  xp->cooling_data.e_frac = xp->cooling_data.HII_frac \
    + 0.25 * xp->cooling_data.HeII_frac \
    + 0.5  * xp->cooling_data.HeIII_frac;
#endif  // MODE >= 1
lhausamm's avatar
lhausamm committed
212
213
214

#if COOLING_GRACKLE_MODE >= 2
  /* primordial chemistry >= 2 */
lhausamm's avatar
lhausamm committed
215
216
217
  xp->cooling_data.HM_frac = zero;
  xp->cooling_data.H2I_frac = zero;
  xp->cooling_data.H2II_frac = zero;
218
#endif  // MODE >= 2
lhausamm's avatar
lhausamm committed
219
220
221

#if COOLING_GRACKLE_MODE >= 3
  /* primordial chemistry >= 3 */
222
223
  xp->cooling_data.DI_frac = grackle_data->DeuteriumToHydrogenRatio
    * grackle_data->HydrogenFractionByMass;
lhausamm's avatar
lhausamm committed
224
225
  xp->cooling_data.DII_frac = zero;
  xp->cooling_data.HDI_frac = zero;
lhausamm's avatar
lhausamm committed
226
#endif  // MODE >= 3
lhausamm's avatar
lhausamm committed
227

228
#if COOLING_GRACKLE_MODE > 0
lhausamm's avatar
lhausamm committed
229
  cooling_compute_equilibrium(p, xp, cooling);
230
#endif
231
232
}

233
234
235
236
237
238
239
240
241
242
243
/**
 * @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;
}

244
245
246
247
248
/**
 * @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
249
__attribute__((always_inline)) INLINE static void cooling_print_backend(
250
251
252
    const struct cooling_function_data* cooling) {

  message("Cooling function is 'Grackle'.");
253
254
  message("Using Grackle    = %i", cooling->chemistry.use_grackle);
  message("Chemical network = %i",
Matthieu Schaller's avatar
Matthieu Schaller committed
255
          cooling->chemistry.primordial_chemistry);
256
257
258
259
260
261
262
263
264
  message("CloudyTable      = %s", cooling->cloudy_table);
  message("Redshift         = %g", cooling->redshift);
  message("UV background    = %d", cooling->with_uv_background);
  message("Metal cooling    = %i", cooling->chemistry.metal_cooling);
  message("Self Shielding   = %i", cooling->self_shielding_method);
  message("Specific Heating Rates   = %i",
	  cooling->provide_specific_heating_rates);
  message("Volumetric Heating Rates = %i",
	  cooling->provide_volumetric_heating_rates);
265
266
  message("Units:");
  message("\tComoving     = %i", cooling->units.comoving_coordinates);
lhausamm's avatar
lhausamm committed
267
268
  message("\tLength       = %g", cooling->units.length_units);
  message("\tDensity      = %g", cooling->units.density_units);
269
  message("\tTime         = %g", cooling->units.time_units);
lhausamm's avatar
lhausamm committed
270
  message("\tScale Factor = %g", cooling->units.a_units);
271
272
}

lhausamm's avatar
lhausamm committed
273
/**
274
275
276
277
278
279
 * @brief copy a single field from the grackle data to a #xpart
 * 
 * @param data The #grackle_field_data
 * @param xp The #xpart
 * @param rho Particle density
 * @param field The field to copy
lhausamm's avatar
lhausamm committed
280
 */
281
282
#define cooling_copy_field_from_grackle(data, xp, rho, field)		\
  xp->cooling_data.field ## _frac = *data.field ## _density / rho;
lhausamm's avatar
lhausamm committed
283

284
285
286
287
288
289
290
291
292
293
294
/**
 * @brief copy a single field from a #xpart to the grackle data
 * 
 * @param data The #grackle_field_data
 * @param xp The #xpart
 * @param rho Particle density
 * @param field The field to copy
 */
#define cooling_copy_field_to_grackle(data, xp, rho, field)		\
  gr_float grackle_ ## field = xp->cooling_data.field ## _frac * rho;	\
  data.field ## _density = &grackle_ ## field;
lhausamm's avatar
lhausamm committed
295
296
297


/**
298
 * @brief copy a #xpart to the grackle data
lhausamm's avatar
lhausamm committed
299
 *
300
301
302
303
304
305
306
 * Warning this function creates some variable, therefore the grackle call
 * should be in a block that still has the variables.
 * 
 * @param data The #grackle_field_data
 * @param p The #part
 * @param xp The #xpart
 * @param rho Particle density
lhausamm's avatar
lhausamm committed
307
 */
308
309
310
311
312
313
314
315
316
#if COOLING_GRACKLE_MODE > 0
#define cooling_copy_to_grackle1(data, p, xp, rho)			\
  cooling_copy_field_to_grackle(data, xp, rho, HI);			\
  cooling_copy_field_to_grackle(data, xp, rho, HII);			\
  cooling_copy_field_to_grackle(data, xp, rho, HeI);			\
  cooling_copy_field_to_grackle(data, xp, rho, HeII);			\
  cooling_copy_field_to_grackle(data, xp, rho, HeIII);			\
  cooling_copy_field_to_grackle(data, xp, rho, e);
#else
317
318
319
320
321
322
323
#define cooling_copy_to_grackle1(data, p, xp, rho)	\
  data.HI_density = NULL;				\
  data.HII_density = NULL;				\
  data.HeI_density = NULL;				\
  data.HeII_density = NULL;				\
  data.HeIII_density = NULL;				\
  data.e_density = NULL;
324
#endif
lhausamm's avatar
lhausamm committed
325
326

/**
327
 * @brief copy a #xpart to the grackle data
lhausamm's avatar
lhausamm committed
328
 *
329
330
331
332
333
334
335
 * Warning this function creates some variable, therefore the grackle call
 * should be in a block that still has the variables.
 * 
 * @param data The #grackle_field_data
 * @param p The #part
 * @param xp The #xpart
 * @param rho Particle density
lhausamm's avatar
lhausamm committed
336
 */
337
338
339
340
341
342
#if COOLING_GRACKLE_MODE > 1
#define cooling_copy_to_grackle2(data, p, xp, rho)			\
  cooling_copy_field_to_grackle(data, xp, rho, HM);			\
  cooling_copy_field_to_grackle(data, xp, rho, H2I);			\
  cooling_copy_field_to_grackle(data, xp, rho, H2II);
#else
343
344
345
346
#define cooling_copy_to_grackle2(data, p, xp, rho)	\
  data.HM_density = NULL;				\
  data.H2I_density = NULL;				\
  data.H2II_density = NULL;
347
#endif
lhausamm's avatar
lhausamm committed
348
349

/**
350
 * @brief copy a #xpart to the grackle data
lhausamm's avatar
lhausamm committed
351
 *
352
353
354
355
356
357
358
 * Warning this function creates some variable, therefore the grackle call
 * should be in a block that still has the variables.
 * 
 * @param data The #grackle_field_data
 * @param p The #part
 * @param xp The #xpart
 * @param rho Particle density
lhausamm's avatar
lhausamm committed
359
 */
360
361
362
363
364
365
#if COOLING_GRACKLE_MODE > 2
#define cooling_copy_to_grackle3(data, p, xp, rho)			\
  cooling_copy_field_to_grackle(data, xp, rho, DI);			\
  cooling_copy_field_to_grackle(data, xp, rho, DII);			\
  cooling_copy_field_to_grackle(data, xp, rho, HDI);
#else
366
367
368
369
#define cooling_copy_to_grackle3(data, p, xp, rho)	\
  data.DI_density = NULL;				\
  data.DII_density = NULL;				\
  data.HDI_density = NULL;
370
#endif
lhausamm's avatar
lhausamm committed
371

372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
/**
 * @brief copy the grackle data to a #xpart
 * 
 * @param data The #grackle_field_data
 * @param p The #part
 * @param xp The #xpart
 * @param rho Particle density
 */
#if COOLING_GRACKLE_MODE > 0
#define cooling_copy_from_grackle1(data, p, xp, rho)			\
  cooling_copy_field_from_grackle(data, xp, rho, HI);			\
  cooling_copy_field_from_grackle(data, xp, rho, HII);			\
  cooling_copy_field_from_grackle(data, xp, rho, HeI);			\
  cooling_copy_field_from_grackle(data, xp, rho, HeII);			\
  cooling_copy_field_from_grackle(data, xp, rho, HeIII);		\
  cooling_copy_field_from_grackle(data, xp, rho, e);
#else
#define cooling_copy_from_grackle1(data, p, xp, rho)
#endif
lhausamm's avatar
lhausamm committed
391

392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
/**
 * @brief copy the grackle data to a #xpart
 * 
 * @param data The #grackle_field_data
 * @param p The #part
 * @param xp The #xpart
 * @param rho Particle density
 */
#if COOLING_GRACKLE_MODE > 1
#define cooling_copy_from_grackle2(data, p, xp, rho)			\
  cooling_copy_field_from_grackle(data, xp, rho, HM);			\
  cooling_copy_field_from_grackle(data, xp, rho, H2I);			\
  cooling_copy_field_from_grackle(data, xp, rho, H2II);
#else
#define cooling_copy_from_grackle2(data, p, xp, rho)
#endif
lhausamm's avatar
lhausamm committed
408

409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
/**
 * @brief copy the grackle data to a #xpart
 * 
 * @param data The #grackle_field_data
 * @param p The #part
 * @param xp The #xpart
 * @param rho Particle density
 */
#if COOLING_GRACKLE_MODE > 2
#define cooling_copy_from_grackle3(data, p, xp, rho)			\
  cooling_copy_field_from_grackle(data, xp, rho, DI);			\
  cooling_copy_field_from_grackle(data, xp, rho, DII);			\
  cooling_copy_field_from_grackle(data, xp, rho, HDI);
#else
#define cooling_copy_from_grackle3(data, p, xp, rho)
#endif
lhausamm's avatar
lhausamm committed
425
426


427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
/**
 * @brief copy a #xpart to the grackle data
 *
 * Warning this function creates some variable, therefore the grackle call
 * should be in a block that still has the variables.
 * 
 * @param data The #grackle_field_data
 * @param p The #part
 * @param xp The #xpart
 * @param rho Particle density
 */
#define cooling_copy_to_grackle(data, p, xp, rho)			\
  cooling_copy_to_grackle1(data, p, xp, rho);				\
  cooling_copy_to_grackle2(data, p, xp, rho);				\
  cooling_copy_to_grackle3(data, p, xp, rho);				\
442
443
444
445
446
447
448
  data.volumetric_heating_rate = NULL;					\
  data.specific_heating_rate = NULL;					\
  data.RT_heating_rate = NULL;						\
  data.RT_HI_ionization_rate = NULL;					\
  data.RT_HeI_ionization_rate = NULL;					\
  data.RT_HeII_ionization_rate = NULL;					\
  data.RT_H2_dissociation_rate = NULL;					\
449
450
  gr_float metal_density = chemistry_metal_mass_fraction(p, xp) * rho;	\
  data.metal_density = &metal_density;
lhausamm's avatar
lhausamm committed
451

452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
/**
 * @brief copy a #xpart to the grackle data
 *
 * Warning this function creates some variable, therefore the grackle call
 * should be in a block that still has the variables.
 * 
 * @param data The #grackle_field_data
 * @param p The #part
 * @param xp The #xpart
 * @param rho Particle density
 */
#define cooling_copy_from_grackle(data, p, xp, rho)			\
  cooling_copy_from_grackle1(data, p, xp, rho);				\
  cooling_copy_from_grackle2(data, p, xp, rho);				\
  cooling_copy_from_grackle3(data, p, xp, rho);
lhausamm's avatar
lhausamm committed
467

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

lhausamm's avatar
lhausamm committed
487
  /* set current time */
lhausamm's avatar
lhausamm committed
488
  code_units units = cooling->units;
lhausamm's avatar
lhausamm committed
489
490
491
  if (cooling->redshift == -1)
    error("TODO time dependant redshift");
  else
lhausamm's avatar
lhausamm committed
492
    units.a_value = 1. / (1. + cooling->redshift);
493

lhausamm's avatar
lhausamm committed
494
495
  /* initialize data */
  grackle_field_data data;
lhausamm's avatar
lhausamm committed
496

lhausamm's avatar
lhausamm committed
497
498
499
500
501
  /* 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
502

503
  data.grid_dx = 0.;
lhausamm's avatar
lhausamm committed
504
505
506
507
  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
508

lhausamm's avatar
lhausamm committed
509
  /* general particle data */
510
511
  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
512
  gr_float energy = energy_before;
513

lhausamm's avatar
lhausamm committed
514
  /* initialize density */
lhausamm's avatar
lhausamm committed
515
  data.density = &density;
lhausamm's avatar
lhausamm committed
516
517

  /* initialize energy */
lhausamm's avatar
lhausamm committed
518
  data.internal_energy = &energy;
lhausamm's avatar
lhausamm committed
519
520
521
522
523

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

525
526
  /* copy to grackle structure */
  cooling_copy_to_grackle(data, p, xp, density);
lhausamm's avatar
lhausamm committed
527
528

  /* solve chemistry with table */
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551

  
  chemistry_data chemistry_grackle = cooling->chemistry;
  chemistry_data_storage my_rates = grackle_rates;
  _solve_chemistry(&chemistry_grackle,
		   &my_rates,
		   &units, dt, data.grid_dx,
		   data.grid_rank, data.grid_dimension,
		   data.grid_start, data.grid_end,
		   data.density, data.internal_energy,
		   data.x_velocity, data.y_velocity, data.z_velocity,
		   data.HI_density, data.HII_density, data.HM_density,
		   data.HeI_density, data.HeII_density, data.HeIII_density,
		   data.H2I_density, data.H2II_density,
		   data.DI_density, data.DII_density, data.HDI_density,
		   data.e_density, data.metal_density,
		   data.volumetric_heating_rate, data.specific_heating_rate,
		   data.RT_heating_rate, data.RT_HI_ionization_rate, data.RT_HeI_ionization_rate,
		   data.RT_HeII_ionization_rate, data.RT_H2_dissociation_rate,
		   NULL); 
  //if (solve_chemistry(&units, &data, dt) == 0) {
  //  error("Error in solve_chemistry.");
  //}
lhausamm's avatar
lhausamm committed
552

lhausamm's avatar
lhausamm committed
553
  /* copy from grackle data to particle */
554
  cooling_copy_from_grackle(data, p, xp, density);
lhausamm's avatar
lhausamm committed
555

lhausamm's avatar
lhausamm committed
556
  /* compute rate */
lhausamm's avatar
lhausamm committed
557
  return (energy - energy_before) / dt;
558
559
}

560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
/**
 * @brief Compute the cooling time
 *
 * @param cooling The #cooling_function_data used in the run.
 * @param p Pointer to the particle data.
 * @param xp Pointer to the particle extra data
 *
 * @return cooling time
 */
__attribute__((always_inline)) INLINE static gr_float cooling_time(
    const struct cooling_function_data* restrict cooling,
    const struct part* restrict p, struct xpart* restrict xp) {

  /* set current time */
  code_units units = cooling->units;
  if (cooling->redshift == -1)
    error("TODO time dependant redshift");
  else
    units.a_value = 1. / (1. + cooling->redshift);

  /* initialize data */
  grackle_field_data data;

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

  /* general particle data */
  const gr_float energy_before = hydro_get_internal_energy(p);
  gr_float density = hydro_get_density(p);
  gr_float energy = energy_before;

  /* initialize density */
  data.density = &density;

  /* initialize energy */
  data.internal_energy = &energy;

  /* grackle 3.0 doc: "Currently not used" */
  data.x_velocity = NULL;
  data.y_velocity = NULL;
  data.z_velocity = NULL;

  /* copy data from particle to grackle data */
611
  cooling_copy_to_grackle(data, p, xp, density);
612
613
614
615
616
617
618
619

  /* Compute cooling time */
  gr_float cooling_time;
  if (calculate_cooling_time(&units, &data, &cooling_time) == 0) {
    error("Error in calculate_cooling_time.");
  }

  /* copy from grackle data to particle */
620
  cooling_copy_from_grackle(data, p, xp, density);
621
622
623
624
625

  /* compute rate */
  return cooling_time;
}

626
627
628
629
630
/**
 * @brief Apply the cooling function to a particle.
 *
 * @param phys_const The physical constants in internal units.
 * @param us The internal system of units.
631
 * @param cosmo The current cosmological model.
632
633
634
635
636
637
 * @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
638
    const struct unit_system* restrict us,
639
    const struct cosmology* restrict cosmo,
640
    const struct cooling_function_data* restrict cooling,
lhausamm's avatar
lhausamm committed
641
    struct part* restrict p, struct xpart* restrict xp, double dt) {
642

Matthieu Schaller's avatar
Matthieu Schaller committed
643
  if (dt == 0.) return;
644
645
646

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

648
  /* compute cooling rate */
649
  const float du_dt = cooling_rate(phys_const, us, cosmo, cooling, p, dt);
650
651

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

654
  /* Update the internal energy */
655
  hydro_set_internal_energy_dt(p, hydro_du_dt + du_dt);
656
657
658
659
660
661
662
663
664
}

/**
 * @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.
665
 * @param cosmo The current cosmological model.
666
667
668
669
670
671
 * @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,
672
    const struct cosmology* restrict cosmo,
lhausamm's avatar
lhausamm committed
673
    const struct unit_system* restrict us, const struct part* restrict p) {
674
675
676
677

  return FLT_MAX;
}

678
679
680
681
682
683
684
685
686
687
688
689
690


/**
 * @brief Initialises the cooling unit system.
 *
 * @param us The current internal system of units.
 * @param cooling The cooling properties to initialize
 */
__attribute__((always_inline)) INLINE static void cooling_init_units(
    const struct unit_system* us,
    struct cooling_function_data* cooling) {

  /* These are conversions from code units to cgs. */
691
692
693

  /* first cosmo */
  cooling->units.a_units = 1.0;  // units for the expansion factor (1/1+zi)
694
  cooling->units.a_value = 1.0;
695
696
697

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

  /* then units */
Matthieu Schaller's avatar
Matthieu Schaller committed
701
702
  cooling->units.density_units =
      us->UnitMass_in_cgs / pow(us->UnitLength_in_cgs, 3);
703
704
  cooling->units.length_units = us->UnitLength_in_cgs;
  cooling->units.time_units = us->UnitTime_in_cgs;
Matthieu Schaller's avatar
Matthieu Schaller committed
705
706
707
  cooling->units.velocity_units = cooling->units.a_units *
                                  cooling->units.length_units /
                                  cooling->units.time_units;
708
709
710
711
712
713
714
715
716
717
718
719
720
721
}

/**
 * @brief Initialises Grackle.
 *
 * @param cooling The cooling properties to initialize
 */
__attribute__((always_inline)) INLINE static void cooling_init_grackle(
    struct cooling_function_data* cooling) {
  
#ifdef SWIFT_DEBUG_CHECKS
  /* enable verbose for grackle */
  grackle_verbose = 1;
#endif
lhausamm's avatar
lhausamm committed
722

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

725
  /* Create a chemistry object for parameters and rate data. */
lhausamm's avatar
lhausamm committed
726
  if (set_default_chemistry_parameters(chemistry) == 0) {
727
728
    error("Error in set_default_chemistry_parameters.");
  }
lhausamm's avatar
lhausamm committed
729

730
  // Set parameter values for chemistry.
lhausamm's avatar
lhausamm committed
731
732
  chemistry->use_grackle = 1;
  chemistry->with_radiative_cooling = 1;
lhausamm's avatar
lhausamm committed
733

734
735
  /* molecular network with H, He, D
   From Cloudy table */
736
  chemistry->primordial_chemistry = COOLING_GRACKLE_MODE;
737
738
  chemistry->metal_cooling = cooling->with_metal_cooling;
  chemistry->UVbackground = cooling->with_uv_background;
lhausamm's avatar
lhausamm committed
739
  chemistry->grackle_data_file = cooling->cloudy_table;
lhausamm's avatar
lhausamm committed
740

741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
  /* radiative transfer */
  chemistry->use_radiative_transfer =
    cooling->provide_specific_heating_rates ||
    cooling->provide_volumetric_heating_rates;
  chemistry->use_volumetric_heating_rate =
    cooling->provide_volumetric_heating_rates;
  chemistry->use_specific_heating_rate =
    cooling->provide_specific_heating_rates;

  if (cooling->provide_specific_heating_rates &&
      cooling->provide_volumetric_heating_rates)
    message("WARNING: You should specified either the specific or the volumetric heating rates, not both");

  /* self shielding */
  chemistry->self_shielding_method =
    cooling->self_shielding_method;
lhausamm's avatar
lhausamm committed
757

lhausamm's avatar
lhausamm committed
758
759
  /* Initialize the chemistry object. */
  if (initialize_chemistry_data(&cooling->units) == 0) {
760
761
    error("Error in initialize_chemistry_data.");
  }
lhausamm's avatar
lhausamm committed
762

763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
}
  
/**
 * @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
 */
__attribute__((always_inline)) INLINE static void cooling_init_backend(
    const struct swift_params* parameter_file, const struct unit_system* us,
    const struct phys_const* phys_const,
    struct cooling_function_data* cooling) {

  if (GRACKLE_NPART != 1)
    error("Grackle with multiple particles not implemented");

  /* read parameters */
782
  cooling_read_parameters(parameter_file, cooling);
783
784
785
786
787
788

  /* Set up the units system. */
  cooling_init_units(us, cooling);

  cooling_init_grackle(cooling);

789
790
791
}

#endif /* SWIFT_COOLING_GRACKLE_H */