logger_particle.c 6.85 KB
Newer Older
Loic Hausammann's avatar
Loic Hausammann committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
 *
 * 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
#include "logger_particle.h"
#include "logger_header.h"
#include "logger_io.h"
Loic Hausammann's avatar
Loic Hausammann committed
22
#include "logger_reader.h"
23
#include "logger_time.h"
Loic Hausammann's avatar
Loic Hausammann committed
24
#include "logger_tools.h"
25

Loic Hausammann's avatar
format    
Loic Hausammann committed
26
#include <stdint.h>
27
28
29
#include <stdio.h>
#include <stdlib.h>

30
/**
Loic Hausammann's avatar
Loic Hausammann committed
31
 * @brief Print the properties of a logger_particle.
32
 *
Loic Hausammann's avatar
Loic Hausammann committed
33
 * @param p The #logger_particle to print
34
 */
Loic Hausammann's avatar
Loic Hausammann committed
35
36
37
38
39
40
41
42
43
44
void logger_particle_print(const struct logger_particle *p) {
  message("ID:            %lu\n", p->id);
  message("Mass:          %g\n", p->mass);
  message("Time:          %g\n", p->time);
  message("Cutoff Radius: %g\n", p->h);
  message("Positions:     (%g, %g, %g)\n", p->pos[0], p->pos[1], p->pos[2]);
  message("Velocities:    (%g, %g, %g)\n", p->vel[0], p->vel[1], p->vel[2]);
  message("Accelerations: (%g, %g, %g)\n", p->acc[0], p->acc[1], p->acc[2]);
  message("Entropy:       %g\n", p->entropy);
  message("Density:       %g\n", p->density);
45
46
}

47
/**
Loic Hausammann's avatar
Loic Hausammann committed
48
 * @brief Initialize a logger_particle.
49
 *
Loic Hausammann's avatar
Loic Hausammann committed
50
 * @param part The #logger_particle to initialize.
51
 */
Loic Hausammann's avatar
Loic Hausammann committed
52
void logger_particle_init(struct logger_particle *part) {
53
54
55
56
57
58
59
60
61
62
63
64
65
  for (size_t k = 0; k < DIM; k++) {
    part->pos[k] = 0;
    part->vel[k] = 0;
    part->acc[k] = 0;
  }

  part->entropy = -1;
  part->density = -1;
  part->h = -1;
  part->mass = -1;
  part->id = SIZE_MAX;
}

66
/**
Loic Hausammann's avatar
Loic Hausammann committed
67
 * @brief Read a single field for a particle.
68
 *
Loic Hausammann's avatar
Loic Hausammann committed
69
 * @param part The #logger_particle to update.
Loic Hausammann's avatar
Loic Hausammann committed
70
 * @param data Pointer to the data to read.
Loic Hausammann's avatar
Loic Hausammann committed
71
72
73
 * @param offset position to read.
 * @param field field to read.
 * @param size number of bits to read.
74
 *
Loic Hausammann's avatar
Loic Hausammann committed
75
 * @return position after the data read.
76
 */
Loic Hausammann's avatar
Loic Hausammann committed
77
78
79
size_t logger_particle_read_field(struct logger_particle *part, void *map,
				  size_t offset, const char *field,
				  const size_t size) {
80
81
  void *p = NULL;

Loic Hausammann's avatar
Loic Hausammann committed
82
  if (strcmp("positions", field) == 0) {
83
    p = &part->pos;
Loic Hausammann's avatar
Loic Hausammann committed
84
  } else if (strcmp("velocities", field) == 0) {
85
    p = &part->vel;
Loic Hausammann's avatar
Loic Hausammann committed
86
  } else if (strcmp("accelerations", field) == 0) {
87
88
89
    p = &part->acc;
  } else if (strcmp("entropy", field) == 0) {
    p = &part->entropy;
Loic Hausammann's avatar
Loic Hausammann committed
90
  } else if (strcmp("smoothing length", field) == 0) {
91
92
93
94
95
96
    p = &part->h;
  } else if (strcmp("density", field) == 0) {
    p = &part->density;
  } else if (strcmp("consts", field) == 0) {
    p = malloc(size);
  } else {
97
    error("Type %s not defined", field);
98
99
  }

Loic Hausammann's avatar
Loic Hausammann committed
100
  offset = io_read_data(map, size, p, offset);
101
102
103
104
105
106
107
108
109
110

  if (strcmp("consts", field) == 0) {
    part->mass = 0;
    part->id = 0;
    memcpy(&part->mass, p, sizeof(float));
    p += sizeof(float);
    memcpy(&part->id, p, sizeof(size_t));
    p -= sizeof(float);
    free(p);
  }
Loic Hausammann's avatar
Loic Hausammann committed
111
112

  return offset;
113
114
}

115
/**
Loic Hausammann's avatar
Loic Hausammann committed
116
 * @brief Read a particle in the dump file.
117
 *
Loic Hausammann's avatar
Loic Hausammann committed
118
 * @param reader The #logger_reader.
Loic Hausammann's avatar
Loic Hausammann committed
119
 * @param part The #logger_particle to update.
Loic Hausammann's avatar
Loic Hausammann committed
120
 * @param offset offset of the chunk to read.
Loic Hausammann's avatar
Loic Hausammann committed
121
 * @param time time to interpolate (if #logger_reader_type is an interpolating
Loic Hausammann's avatar
Loic Hausammann committed
122
123
 * one).
 * @param reader_type #logger_reader_type.
124
 *
Loic Hausammann's avatar
Loic Hausammann committed
125
 * @return position after the record.
126
 */
Loic Hausammann's avatar
Loic Hausammann committed
127
128
129
130
131
132
133
134
135
size_t logger_particle_read(struct logger_particle *part, const struct logger_reader *reader,
			    size_t offset, const double time,
			    const enum logger_reader_type reader_type) {

  const struct header *h = &reader->dump.header;
  void *map = reader->dump.dump.map;

  const struct time_array *times = &reader->dump.times;

136
137
138
  size_t mask = 0;
  size_t h_offset = 0;

Loic Hausammann's avatar
Loic Hausammann committed
139
  logger_particle_init(part);
140

Loic Hausammann's avatar
Loic Hausammann committed
141
  offset = io_read_mask(h, map, offset, &mask, &h_offset);
Loic Hausammann's avatar
Loic Hausammann committed
142

143
  if (mask != 127) error("Unexpected mask: %lu", mask);
144

Loic Hausammann's avatar
Loic Hausammann committed
145
146
  for (size_t i = 0; i < h->number_mask; i++) {
    if (mask & h->masks[i].mask) {
Loic Hausammann's avatar
Loic Hausammann committed
147
148
      offset = logger_particle_read_field(part, map, offset, h->masks[i].name,
					  h->masks[i].size);
149
150
151
    }
  }

Loic Hausammann's avatar
Loic Hausammann committed
152
153
154
155
  if (times) {
    /* move offset by 1 in order to be in the required chunk */
    part->time = time_array_get_time(times, offset - 1);
  }
156
157
158
159
  else
    part->time = -1;

  /* end of const case */
Loic Hausammann's avatar
Loic Hausammann committed
160
161
  if (reader_type == logger_reader_const)
    return offset;
162
163

  /* read next particle */
Loic Hausammann's avatar
Loic Hausammann committed
164
  struct logger_particle part_next;
165

Loic Hausammann's avatar
Loic Hausammann committed
166
  if (!header_is_forward(h)) {
Loic Hausammann's avatar
Loic Hausammann committed
167
168
    error("Cannot read a particle with non forward offsets.");
  }
169

Loic Hausammann's avatar
Loic Hausammann committed
170
171
172
  if (h_offset == 0)
    return offset;

173
  /* get absolute offset of next particle */
Loic Hausammann's avatar
Loic Hausammann committed
174
  h_offset += offset - header_get_mask_size(h, mask) - LOGGER_MASK_SIZE -
Loic Hausammann's avatar
format    
Loic Hausammann committed
175
              LOGGER_OFFSET_SIZE;
176
177
178
179

  part_next.time = time_array_get_time(times, h_offset);

  /* previous part exists */
Loic Hausammann's avatar
Loic Hausammann committed
180
181
  h_offset = logger_particle_read(&part_next, reader, h_offset, part_next.time,
				  logger_reader_const);
182

Loic Hausammann's avatar
Loic Hausammann committed
183
  logger_particle_interpolate(part, &part_next, time);
Loic Hausammann's avatar
Loic Hausammann committed
184
185

  return offset;
186
187
}

188
189
190
/**
 * @brief interpolate two particles at a given time
 *
Loic Hausammann's avatar
Loic Hausammann committed
191
 * @param part_curr #logger_particle In: current particle (before time), Out:
192
 * interpolated particle
Loic Hausammann's avatar
Loic Hausammann committed
193
 * @param part_next #logger_particle next particle (after time)
194
195
196
 * @param time interpolation time
 *
 */
Loic Hausammann's avatar
Loic Hausammann committed
197
198
199
void logger_particle_interpolate(struct logger_particle *part_curr,
                                 const struct logger_particle *part_next,
                                 const double time) {
200

201
202
  if (!part_curr) error("part_curr is NULL");
  if (!part_next) error("part_next is NULL");
203
204
205

#ifdef SWIFT_DEBUG_CHECKS
  if (part_next->time <= part_curr->time)
206
    error("Wrong particle order (next before current)");
207
  if ((time < part_curr->time) || (part_next->time < time))
Loic Hausammann's avatar
Loic Hausammann committed
208
209
210
211
    error(
        "Interpolating, not extrapolating (particle time: %f, "
        "interpolating time: %f, next particle time: %f)",
        part_curr->time, time, part_next->time);
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
#endif

  double scaling = part_next->time - part_curr->time;

  scaling = (time - part_curr->time) / scaling;

  double tmp;
  float ftmp;

  /* interpolate vectors */
  for (size_t i = 0; i < DIM; i++) {
    tmp = (part_next->pos[i] - part_curr->pos[i]);
    part_curr->pos[i] += tmp * scaling;

    ftmp = (part_next->vel[i] - part_curr->vel[i]);
    part_curr->vel[i] += ftmp * scaling;

    ftmp = (part_next->acc[i] - part_curr->acc[i]);
    part_curr->acc[i] += ftmp * scaling;
  }

  /* interpolate scalars */
  ftmp = (part_next->entropy - part_curr->entropy);
  part_curr->entropy += ftmp * scaling;

  /* set time */
  part_curr->time = time;
}