hydro_io.h 7.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/*******************************************************************************
 * This file is part of SWIFT.
 * Coypright (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/>.
 *
 ******************************************************************************/
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#ifndef SWIFT_MINIMAL_HYDRO_IO_H
#define SWIFT_MINIMAL_HYDRO_IO_H

/**
 * @file Minimal/hydro_io.h
 * @brief Minimal conservative implementation of SPH (i/o routines)
 *
 * The thermal variable is the internal energy (u). Simple constant
 * viscosity term without switches is implemented. No thermal conduction
 * term is implemented.
 *
 * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
 * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
 * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
 * pp. 759-794.
 */
35

36
37
#include "adiabatic_index.h"
#include "hydro.h"
38
39
40
#include "io_properties.h"
#include "kernel_hydro.h"

41
/**
42
 * @brief Specifies which particle fields to read from a dataset
43
 *
44
45
46
 * @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.
47
 */
48
49
50
void hydro_read_particles(struct part* parts, struct io_props* list,
                          int* num_fields) {

51
52
  *num_fields = 8;

53
54
55
56
57
58
59
60
61
62
  /* 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,
                                parts, mass);
  list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
                                UNIT_CONV_LENGTH, parts, h);
  list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
63
                                UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
64
65
66
67
68
69
  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, rho);
70
71
}

72
73
void convert_S(const struct engine* e, const struct part* p,
               const struct xpart* xp, float* ret) {
74

75
  ret[0] = hydro_get_comoving_entropy(p);
76
77
}

78
79
void convert_P(const struct engine* e, const struct part* p,
               const struct xpart* xp, float* ret) {
80

81
  ret[0] = hydro_get_comoving_pressure(p);
82
83
84
}

void convert_part_pos(const struct engine* e, const struct part* p,
85
                      const struct xpart* xp, double* ret) {
86
87
88
89
90
91
92
93
94
95

  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];
  }
96
97
}

98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
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;
}

132
/**
133
 * @brief Specifies which particle fields to write to a dataset
134
 *
135
 * @param parts The particle array.
136
 * @param xparts The extended particle array.
137
138
 * @param list The list of i/o properties to write.
 * @param num_fields The number of i/o fields to write.
139
 */
140
141
void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
                           struct io_props* list, int* num_fields) {
142

143
  *num_fields = 9;
144

145
  /* List what we want to write */
146
147
148
149
150
  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);
151
152
153
154
  list[2] =
      io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
  list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
                                 parts, h);
155
  list[4] = io_make_output_field("InternalEnergy", FLOAT, 1,
156
157
158
                                 UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
  list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
                                 UNIT_CONV_NO_UNITS, parts, id);
159
  list[6] =
160
      io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
161
162
163
  list[7] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
                                              UNIT_CONV_ENTROPY_PER_UNIT_MASS,
                                              parts, xparts, convert_S);
164
  list[8] = io_make_output_field_convert_part(
165
      "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, xparts, convert_P);
166
}
167
168
169
170
171

/**
 * @brief Writes the current model of SPH to the file
 * @param h_grpsph The HDF5 group in which to write
 */
172
void hydro_write_flavour(hid_t h_grpsph) {
173
174
175

  /* Viscosity and thermal conduction */
  /* Nothing in this minimal model... */
176
177
178
  io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
  io_write_attribute_s(h_grpsph, "Viscosity Model",
                       "Minimal treatment as in Monaghan (1992)");
179
180

  /* Time integration properties */
181
182
  io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
                       const_max_u_change);
183
}
184
185
186
187
188
189
190

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

#endif /* SWIFT_MINIMAL_HYDRO_IO_H */