hydro_io.h 7.67 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.
 * 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/>.
 *
 ******************************************************************************/
#ifndef SWIFT_PRESSURE_ENTROPY_HYDRO_IO_H
#define SWIFT_PRESSURE_ENTROPY_HYDRO_IO_H

/**
 * @file PressureEntropy/hydro_io.h
24
 * @brief Pressure-Entropy implementation of SPH (i/o routines)
25
26
27
28
 *
 * The thermal variable is the entropy (S) and the entropy is smoothed over
 * contact discontinuities to prevent spurious surface tension.
 *
29
30
 * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
 * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
 */

#include "adiabatic_index.h"
#include "hydro.h"
#include "io_properties.h"
#include "kernel_hydro.h"

/**
 * @brief Specifies which particle fields to read from a dataset
 *
 * @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.
 */
void hydro_read_particles(struct part* parts, struct io_props* list,
                          int* num_fields) {

48
49
  *num_fields = 8;

50
51
52
53
54
55
56
57
58
  /* 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);
59
60
61
  list[4] =
      io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
                          UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, entropy);
62
63
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
void convert_u(const struct engine* e, const struct part* p,
               const struct xpart* xp, float* ret) {
72

73
  ret[0] = hydro_get_comoving_internal_energy(p);
74
75
}

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

79
  ret[0] = hydro_get_comoving_pressure(p);
80
81
82
}

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

  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];
  }
94
95
}

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
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;
}

130
131
132
133
134
135
136
/**
 * @brief Specifies which particle fields to write to a dataset
 *
 * @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.
 */
137
138
void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
                           struct io_props* list, int* num_fields) {
139

140
  *num_fields = 10;
141

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

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

  /* Viscosity and thermal conduction */
  /* Nothing in this minimal model... */
175
176
  io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
  io_write_attribute_s(
177
178
      h_grpsph, "Viscosity Model",
      "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
179
180
  io_write_attribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
  io_write_attribute_f(h_grpsph, "Viscosity beta", 3.f);
181
182

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

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

#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_IO_H */