hydro_properties.c 12.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
/*******************************************************************************
 * 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/>.
 *
 ******************************************************************************/

/* This object's header. */
#include "hydro_properties.h"

/* Standard headers */
24
#include <float.h>
25
26
27
#include <math.h>

/* Local headers. */
28
#include "adiabatic_index.h"
29
#include "common_io.h"
30
#include "dimension.h"
31
#include "equation_of_state.h"
32
33
#include "error.h"
#include "hydro.h"
34
#include "kernel_hydro.h"
35

36
#define hydro_props_default_max_iterations 30
37
#define hydro_props_default_volume_change 1.4f
38
#define hydro_props_default_h_max FLT_MAX
39
#define hydro_props_default_h_tolerance 1e-4
40
#define hydro_props_default_init_temp 0.f
41
#define hydro_props_default_min_temp 0.f
42
#define hydro_props_default_H_ionization_temperature 1e4
Josh Borrow's avatar
Josh Borrow committed
43
44
45
46
47
48
49
#define hydro_props_default_viscosity_alpha 0.8f
#define hydro_props_default_viscosity_alpha_min \
  0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
#define hydro_props_default_viscosity_alpha_max \
  2.0f /* Values taken from (Price,2004), not used in legacy gadget mode */
#define hydro_props_default_viscosity_length \
  0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
50

51
52
53
54
55
56
57
58
/**
 * @brief Initialize the global properties of the hydro scheme.
 *
 * @param p The #hydro_props.
 * @param phys_const The physical constants in the internal unit system.
 * @param us The internal unit system.
 * @param params The parsed parameters.
 */
59
void hydro_props_init(struct hydro_props *p,
60
61
                      const struct phys_const *phys_const,
                      const struct unit_system *us,
lhausamm's avatar
lhausamm committed
62
                      struct swift_params *params) {
63
64
65

  /* Kernel properties */
  p->eta_neighbours = parser_get_param_float(params, "SPH:resolution_eta");
66
67
68
69
70
71

  /* Tolerance for the smoothing length Newton-Raphson scheme */
  p->h_tolerance = parser_get_opt_param_float(params, "SPH:h_tolerance",
                                              hydro_props_default_h_tolerance);

  /* Get derived properties */
72
  p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
73
74
75
76
  const float delta_eta = p->eta_neighbours * (1.f + p->h_tolerance);
  p->delta_neighbours =
      (pow_dimension(delta_eta) - pow_dimension(p->eta_neighbours)) *
      kernel_norm;
77

78
79
80
81
#ifdef SHADOWFAX_SPH
  /* change the meaning of target_neighbours and delta_neighbours */
  p->target_neighbours = 1.0f;
  p->delta_neighbours = 0.0f;
82
  p->eta_neighbours = 1.0f;
83
84
#endif

85
86
87
88
89
  /* Maximal smoothing length */
  p->h_max = parser_get_opt_param_float(params, "SPH:h_max",
                                        hydro_props_default_h_max);

  /* Number of iterations to converge h */
90
91
  p->max_smoothing_iterations = parser_get_opt_param_int(
      params, "SPH:max_ghost_iterations", hydro_props_default_max_iterations);
92
93
94

  /* Time integration properties */
  p->CFL_condition = parser_get_param_float(params, "SPH:CFL_condition");
95
96
  const float max_volume_change = parser_get_opt_param_float(
      params, "SPH:max_volume_change", hydro_props_default_volume_change);
97
  p->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
98
99
100
101
102

  /* Initial temperature */
  p->initial_temperature = parser_get_opt_param_float(
      params, "SPH:initial_temperature", hydro_props_default_init_temp);

103
  /* Minimal temperature */
104
105
106
  p->minimal_temperature = parser_get_opt_param_float(
      params, "SPH:minimal_temperature", hydro_props_default_min_temp);

107
108
  if ((p->initial_temperature != 0.) &&
      (p->initial_temperature < p->minimal_temperature))
109
110
    error("Initial temperature lower than minimal allowed temperature!");

111
112
113
114
115
  /* Neutral to ionized Hydrogen transition temperature */
  p->hydrogen_ionization_temperature =
      parser_get_opt_param_double(params, "SPH:H_ionization_temperature",
                                  hydro_props_default_H_ionization_temperature);

116
  /* Hydrogen mass fraction */
117
118
  const float default_H_fraction =
      1. - phys_const->const_primordial_He_fraction;
119
  p->hydrogen_mass_fraction = parser_get_opt_param_double(
120
      params, "SPH:H_mass_fraction", default_H_fraction);
121

Josh Borrow's avatar
Josh Borrow committed
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
  /* Read the artificial viscosity parameters from the file, if they exist */
  p->viscosity.alpha = parser_get_opt_param_float(
    params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha
  );
  
  p->viscosity.alpha_max = parser_get_opt_param_float(
    params, "SPH:viscosity_alpha_max", hydro_props_default_viscosity_alpha_max
  );

  p->viscosity.alpha_min = parser_get_opt_param_float(
    params, "SPH:viscosity_alpha_min", hydro_props_default_viscosity_alpha_min
  );

  p->viscosity.length = parser_get_opt_param_float(
    params, "SPH:viscosity_length", hydro_props_default_viscosity_length
  );

139
  /* Compute the initial energy (Note the temp. read is in internal units) */
140
  /* u_init = k_B T_init / (mu m_p (gamma - 1)) */
141
142
143
144
  double u_init = phys_const->const_boltzmann_k / phys_const->const_proton_mass;
  u_init *= p->initial_temperature;
  u_init *= hydro_one_over_gamma_minus_one;

145
  /* Correct for hydrogen mass fraction (mu) */
146
  double mean_molecular_weight;
147
  if (p->initial_temperature > p->hydrogen_ionization_temperature)
148
149
150
151
152
    mean_molecular_weight = 4. / (8. - 5. * (1. - p->hydrogen_mass_fraction));
  else
    mean_molecular_weight = 4. / (1. + 3. * p->hydrogen_mass_fraction);

  p->initial_internal_energy = u_init / mean_molecular_weight;
153
154

  /* Compute the minimal energy (Note the temp. read is in internal units) */
155
  /* u_min = k_B T_min / (mu m_p (gamma - 1)) */
156
  double u_min = phys_const->const_boltzmann_k / phys_const->const_proton_mass;
157
  u_min *= p->minimal_temperature;
158
159
  u_min *= hydro_one_over_gamma_minus_one;

160
  /* Correct for hydrogen mass fraction (mu) */
161
  if (p->minimal_temperature > p->hydrogen_ionization_temperature)
162
163
164
165
166
    mean_molecular_weight = 4. / (8. - 5. * (1. - p->hydrogen_mass_fraction));
  else
    mean_molecular_weight = 4. / (1. + 3. * p->hydrogen_mass_fraction);

  p->minimal_internal_energy = u_min / mean_molecular_weight;
167
168
}

169
170
171
172
173
/**
 * @brief Print the global properties of the hydro scheme.
 *
 * @param p The #hydro_props.
 */
174
175
void hydro_props_print(const struct hydro_props *p) {

176
177
  /* Print equation of state first */
  eos_print(&eos);
178

179
  /* Now SPH */
180
181
182
  message("Hydrodynamic scheme: %s in %dD.", SPH_IMPLEMENTATION,
          (int)hydro_dimension);

183
184
185
  message("Hydrodynamic kernel: %s with eta=%f (%.2f neighbours).", kernel_name,
          p->eta_neighbours, p->target_neighbours);

186
  message("Hydrodynamic relative tolerance in h: %.5f (+/- %.4f neighbours).",
187
          p->h_tolerance, p->delta_neighbours);
188

Peter W. Draper's avatar
Peter W. Draper committed
189
  message("Hydrodynamic integration: CFL parameter: %.4f.", p->CFL_condition);
190

Josh Borrow's avatar
Josh Borrow committed
191
192
193
194
195
  message("Artificial viscosity parameters set to alpha: %.3f, max: %.3f, "
          "min: %.3f, length: %.3f.",
          p->viscosity.alpha, p->viscosity.alpha_max,
          p->viscosity.alpha_min, p->viscosity.length);

196
  message(
Matthieu Schaller's avatar
Matthieu Schaller committed
197
198
      "Hydrodynamic integration: Max change of volume: %.2f "
      "(max|dlog(h)/dt|=%f).",
199
      pow_dimension(expf(p->log_max_h_change)), p->log_max_h_change);
200

201
202
203
  if (p->h_max != hydro_props_default_h_max)
    message("Maximal smoothing length allowed: %.4f", p->h_max);

204
205
206
  if (p->max_smoothing_iterations != hydro_props_default_max_iterations)
    message("Maximal iterations in ghost task set to %d (default is %d)",
            p->max_smoothing_iterations, hydro_props_default_max_iterations);
207
208
209

  if (p->initial_temperature != hydro_props_default_init_temp)
    message("Initial gas temperature set to %f", p->initial_temperature);
210
211

  if (p->minimal_temperature != hydro_props_default_min_temp)
212
    message("Minimal gas temperature set to %f", p->minimal_temperature);
213

Josh Borrow's avatar
Josh Borrow committed
214

Peter W. Draper's avatar
Peter W. Draper committed
215
    // Matthieu: Temporary location for this i/o business.
216

217
#ifdef PLANETARY_SPH
218
219
#ifdef PLANETARY_SPH_NO_BALSARA
  message("Planetary SPH: Balsara switch DISABLED");
220
#else
221
222
223
  message("Planetary SPH: Balsara switch ENABLED");
#endif
#endif
224
}
225
226
227
228

#if defined(HAVE_HDF5)
void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {

229
230
  eos_print_snapshot(h_grpsph, &eos);

231
232
233
234
235
236
  io_write_attribute_i(h_grpsph, "Dimension", (int)hydro_dimension);
  io_write_attribute_s(h_grpsph, "Scheme", SPH_IMPLEMENTATION);
  io_write_attribute_s(h_grpsph, "Kernel function", kernel_name);
  io_write_attribute_f(h_grpsph, "Kernel target N_ngb", p->target_neighbours);
  io_write_attribute_f(h_grpsph, "Kernel delta N_ngb", p->delta_neighbours);
  io_write_attribute_f(h_grpsph, "Kernel eta", p->eta_neighbours);
237
  io_write_attribute_f(h_grpsph, "Smoothing length tolerance", p->h_tolerance);
238
  io_write_attribute_f(h_grpsph, "Maximal smoothing length", p->h_max);
239
240
241
242
243
244
245
  io_write_attribute_f(h_grpsph, "CFL parameter", p->CFL_condition);
  io_write_attribute_f(h_grpsph, "Volume log(max(delta h))",
                       p->log_max_h_change);
  io_write_attribute_f(h_grpsph, "Volume max change time-step",
                       pow_dimension(expf(p->log_max_h_change)));
  io_write_attribute_i(h_grpsph, "Max ghost iterations",
                       p->max_smoothing_iterations);
246
  io_write_attribute_f(h_grpsph, "Minimal temperature", p->minimal_temperature);
247
248
249
250
251
  io_write_attribute_f(h_grpsph, "Initial temperature", p->initial_temperature);
  io_write_attribute_f(h_grpsph, "Initial energy per unit mass",
                       p->initial_internal_energy);
  io_write_attribute_f(h_grpsph, "Hydrogen mass fraction",
                       p->hydrogen_mass_fraction);
252
  io_write_attribute_f(h_grpsph, "Hydrogen ionization transition temperature",
253
                       p->hydrogen_ionization_temperature);
Josh Borrow's avatar
Josh Borrow committed
254
255
256
257
258
259
260
  io_write_attribute_f(h_grpsph, "Alpha viscosity", p->viscosity.alpha);
  io_write_attribute_f(h_grpsph, "Alpha viscosity (max)",
                       p->viscosity.alpha_max);
  io_write_attribute_f(h_grpsph, "Alpha viscosity (min)",
                       p->viscosity.alpha_min);
  io_write_attribute_f(h_grpsph, "Viscosity decay length", p->viscosity.length);
  io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
261
262
}
#endif
263

Josh Borrow's avatar
Josh Borrow committed
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
/**
 * @brief Initialises a hydro_props struct with somewhat useful values for
 *        the automated test suite. This is not intended for production use,
 *        but rather to fill for the purposes of mocking.
 * 
 * @param p the struct
 */
void hydro_props_init_no_hydro(struct hydro_props *p) {
  p->eta_neighbours = 1.2348;
  p->h_tolerance = hydro_props_default_h_tolerance;
  p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
  const float delta_eta = p->eta_neighbours * (1.f + p->h_tolerance);
  p->delta_neighbours =
      (pow_dimension(delta_eta) - pow_dimension(p->eta_neighbours)) *
      kernel_norm;
  p->h_max = hydro_props_default_h_max;
  p->max_smoothing_iterations = hydro_props_default_max_iterations;
  p->CFL_condition = 0.1;
  p->log_max_h_change = logf(powf(1.4, hydro_dimension_inv));
  
  /* These values are inconsistent and in a production run would probably lead
     to a crash. Again, this function is intended for mocking use in unit tests
     and is _not_ to be used otherwise! */
  p->minimal_temperature = hydro_props_default_min_temp;
  p->minimal_internal_energy = 0.f;
  p->initial_temperature = hydro_props_default_init_temp;
  p->initial_internal_energy = 0.f;

  p->hydrogen_mass_fraction = 0.755;
  p->hydrogen_ionization_temperature = hydro_props_default_H_ionization_temperature;

  p->viscosity.alpha = hydro_props_default_viscosity_alpha;
  p->viscosity.alpha_max= hydro_props_default_viscosity_alpha_max;
  p->viscosity.alpha_min= hydro_props_default_viscosity_alpha_min;
  p->viscosity.length= hydro_props_default_viscosity_length;
}

301
302
303
304
305
306
307
/**
 * @brief Write a hydro_props struct to the given FILE as a stream of bytes.
 *
 * @param p the struct
 * @param stream the file stream
 */
void hydro_props_struct_dump(const struct hydro_props *p, FILE *stream) {
Peter W. Draper's avatar
Peter W. Draper committed
308
  restart_write_blocks((void *)p, sizeof(struct hydro_props), 1, stream,
309
                       "hydroprops", "hydro props");
310
311
312
313
314
315
316
317
318
319
}

/**
 * @brief Restore a hydro_props struct from the given FILE as a stream of
 * bytes.
 *
 * @param p the struct
 * @param stream the file stream
 */
void hydro_props_struct_restore(const struct hydro_props *p, FILE *stream) {
Peter W. Draper's avatar
Peter W. Draper committed
320
321
  restart_read_blocks((void *)p, sizeof(struct hydro_props), 1, stream, NULL,
                      "hydro props");
322
}