hydro_io.h 8.71 KB
Newer Older
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Coypright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 *
 * 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/>.
 *
 ******************************************************************************/
19
20
#ifndef SWIFT_GIZMO_HYDRO_IO_H
#define SWIFT_GIZMO_HYDRO_IO_H
21

22
#include "adiabatic_index.h"
23
#include "hydro.h"
24
25
#include "hydro_gradients.h"
#include "hydro_slope_limiters.h"
26
#include "io_properties.h"
27
#include "riemann.h"
28

29
30
31
32
33
34
35
/* Set the description of the particle movement. */
#if defined(GIZMO_FIX_PARTICLES)
#define GIZMO_PARTICLE_MOVEMENT "Fixed particles."
#else
#define GIZMO_PARTICLE_MOVEMENT "Particles move with flow velocity."
#endif

36
/**
37
 * @brief Specifies which particle fields to read from a dataset
38
 *
39
40
41
 * @param parts The particle array.
 * @param list The list of i/o properties to read.
 * @param num_fields The number of i/o fields to read.
42
 */
43
44
45
void hydro_read_particles(struct part* parts, struct io_props* list,
                          int* num_fields) {

46
47
  *num_fields = 8;

48
49
50
51
52
53
  /* List what we want to read */
  list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
                                UNIT_CONV_LENGTH, parts, x);
  list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
                                UNIT_CONV_SPEED, parts, v);
  list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
54
                                parts, conserved.mass);
55
56
  list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
                                UNIT_CONV_LENGTH, parts, h);
57
58
59
  list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
                                UNIT_CONV_ENERGY_PER_UNIT_MASS, parts,
                                conserved.energy);
60
61
62
63
64
65
66
67
  list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
                                UNIT_CONV_NO_UNITS, parts, id);
  list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
                                UNIT_CONV_ACCELERATION, parts, a_hydro);
  list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
                                UNIT_CONV_DENSITY, parts, primitives.rho);
}

68
69
70
71
72
/**
 * @brief Get the internal energy of a particle
 *
 * @param e #engine.
 * @param p Particle.
73
 * @param ret (return) Internal energy of the particle
74
 */
75
76
void convert_u(const struct engine* e, const struct part* p,
               const struct xpart* xp, float* ret) {
77

78
  ret[0] = hydro_get_comoving_internal_energy(p);
79
80
}

81
82
83
84
85
/**
 * @brief Get the entropic function of a particle
 *
 * @param e #engine.
 * @param p Particle.
86
 * @param ret (return) Entropic function of the particle
87
 */
88
89
void convert_A(const struct engine* e, const struct part* p,
               const struct xpart* xp, float* ret) {
90
  ret[0] = hydro_get_comoving_entropy(p);
91
92
}

93
94
95
96
97
98
99
/**
 * @brief Get the total energy of a particle
 *
 * @param e #engine.
 * @param p Particle.
 * @return Total energy of the particle
 */
100
101
void convert_Etot(const struct engine* e, const struct part* p,
                  const struct xpart* xp, float* ret) {
102
#ifdef GIZMO_TOTAL_ENERGY
103
  ret[0] = p->conserved.energy;
104
#else
105
106
107
108
109
110
  float momentum2;

  momentum2 = p->conserved.momentum[0] * p->conserved.momentum[0] +
              p->conserved.momentum[1] * p->conserved.momentum[1] +
              p->conserved.momentum[2] * p->conserved.momentum[2];

111
  ret[0] = p->conserved.energy + 0.5f * momentum2 / p->conserved.mass;
112
#endif
113
114
}

115
void convert_part_pos(const struct engine* e, const struct part* p,
116
                      const struct xpart* xp, double* ret) {
117
118
119
120
121
122
123
124
125
126
127
128

  if (e->s->periodic) {
    ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
    ret[1] = box_wrap(p->x[1], 0.0, e->s->dim[1]);
    ret[2] = box_wrap(p->x[2], 0.0, e->s->dim[2]);
  } else {
    ret[0] = p->x[0];
    ret[1] = p->x[1];
    ret[2] = p->x[2];
  }
}

129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
void convert_part_vel(const struct engine* e, const struct part* p,
                      const struct xpart* xp, float* ret) {

  const int with_cosmology = (e->policy & engine_policy_cosmology);
  const struct cosmology* cosmo = e->cosmology;
  const integertime_t ti_current = e->ti_current;
  const double time_base = e->time_base;

  const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
  const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);

  /* Get time-step since the last kick */
  float dt_kick_grav, dt_kick_hydro;
  if (with_cosmology) {
    dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
    dt_kick_grav -=
        cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
    dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
    dt_kick_hydro -=
        cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
  } else {
    dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
    dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
  }

  /* Extrapolate the velocites to the current time */
  hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);

  /* Conversion from internal units to peculiar velocities */
  ret[0] *= cosmo->a2_inv;
  ret[1] *= cosmo->a2_inv;
  ret[2] *= cosmo->a2_inv;
}

163
/**
164
 * @brief Specifies which particle fields to write to a dataset
165
 *
166
167
168
 * @param parts The particle array.
 * @param list The list of i/o properties to write.
 * @param num_fields The number of i/o fields to write.
169
 */
170
171
void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
                           struct io_props* list, int* num_fields) {
172

173
174
  *num_fields = 10;

175
  /* List what we want to write */
176
177
178
179
180
181
  list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
                                              UNIT_CONV_LENGTH, parts, xparts,
                                              convert_part_pos);
  list[1] = io_make_output_field_convert_part(
      "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);

182
183
184
185
186
187
  list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts,
                                 conserved.mass);
  list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
                                 parts, h);
  list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
                                              UNIT_CONV_ENERGY_PER_UNIT_MASS,
188
                                              parts, xparts, convert_u);
189
190
  list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
                                 UNIT_CONV_NO_UNITS, parts, id);
191
  list[6] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
192
                                 primitives.rho);
193
  list[7] = io_make_output_field_convert_part(
194
      "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, xparts, convert_A);
195
196
  list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
                                 parts, primitives.P);
197
  list[9] = io_make_output_field_convert_part(
198
      "TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, xparts, convert_Etot);
199
200
201
202
203
204
}

/**
 * @brief Writes the current model of SPH to the file
 * @param h_grpsph The HDF5 group in which to write
 */
205
void hydro_write_flavour(hid_t h_grpsph) {
206
  /* Gradient information */
207
208
  io_write_attribute_s(h_grpsph, "Gradient reconstruction model",
                       HYDRO_GRADIENT_IMPLEMENTATION);
209
210

  /* Slope limiter information */
211
212
213
214
  io_write_attribute_s(h_grpsph, "Cell wide slope limiter model",
                       HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION);
  io_write_attribute_s(h_grpsph, "Piecewise slope limiter model",
                       HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION);
215
216

  /* Riemann solver information */
217
218
  io_write_attribute_s(h_grpsph, "Riemann solver type",
                       RIEMANN_SOLVER_IMPLEMENTATION);
219
220
221

  /* Particle movement information */
  io_write_attribute_s(h_grpsph, "Particle movement", GIZMO_PARTICLE_MOVEMENT);
222
}
223
224
225
226
227
228
229

/**
 * @brief Are we writing entropy in the internal energy field ?
 *
 * @return 1 if entropy is in 'internal energy', 0 otherwise.
 */
int writeEntropyFlag() { return 0; }
230
231

#endif /* SWIFT_GIZMO_HYDRO_IO_H */