part.c 9.57 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
24
25
26
27
28
/*******************************************************************************
 * 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/>.
 *
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
#endif

/* This object's header. */
29
#include "error.h"
30
#include "multipole.h"
31
32
#include "part.h"

33
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
34
 * @brief Re-link the #gpart%s associated with the list of #part%s.
35
36
37
 *
 * @param parts The list of #part.
 * @param N The number of particles to re-link;
Matthieu Schaller's avatar
Matthieu Schaller committed
38
 * @param offset The offset of #part%s relative to the global parts list.
39
 */
40
41
void part_relink_gparts_to_parts(struct part *parts, size_t N,
                                 ptrdiff_t offset) {
42
  for (size_t k = 0; k < N; k++) {
43
44
45
    if (parts[k].gpart) {
      parts[k].gpart->id_or_neg_offset = -(k + offset);
    }
46
47
48
49
  }
}

/**
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
 * @brief Re-link the #gpart%s associated with the list of #spart%s.
 *
 * @param sparts The list of #spart.
 * @param N The number of s-particles to re-link;
 * @param offset The offset of #spart%s relative to the global sparts list.
 */
void part_relink_gparts_to_sparts(struct spart *sparts, size_t N,
                                  ptrdiff_t offset) {
  for (size_t k = 0; k < N; k++) {
    if (sparts[k].gpart) {
      sparts[k].gpart->id_or_neg_offset = -(k + offset);
    }
  }
}

/**
 * @brief Re-link the #part%s associated with the list of #gpart%s.
67
68
69
 *
 * @param gparts The list of #gpart.
 * @param N The number of particles to re-link;
70
 * @param parts The global #part array in which to find the #gpart offsets.
71
 */
72
73
void part_relink_parts_to_gparts(struct gpart *gparts, size_t N,
                                 struct part *parts) {
74
  for (size_t k = 0; k < N; k++) {
75
    if (gparts[k].type == swift_type_gas) {
76
77
78
79
80
      parts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
    }
  }
}

81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
/**
 * @brief Re-link the #spart%s associated with the list of #gpart%s.
 *
 * @param gparts The list of #gpart.
 * @param N The number of particles to re-link;
 * @param sparts The global #spart array in which to find the #gpart offsets.
 */
void part_relink_sparts_to_gparts(struct gpart *gparts, size_t N,
                                  struct spart *sparts) {
  for (size_t k = 0; k < N; k++) {
    if (gparts[k].type == swift_type_star) {
      sparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
    }
  }
}

97
98
99
100
101
102
103
104
105
106
107
108
109
/**
 * @brief Verifies that the #gpart, #part and #spart are correctly linked
 * together
 * and that the particle poisitions match.
 *
 * This is a debugging function.
 *
 * @param parts The #part array.
 * @param gparts The #gpart array.
 * @param sparts The #spart array.
 * @param nr_parts The number of #part in the array.
 * @param nr_gparts The number of #gpart in the array.
 * @param nr_sparts The number of #spart in the array.
110
 * @param verbose Do we report verbosely in case of success ?
111
112
113
 */
void part_verify_links(struct part *parts, struct gpart *gparts,
                       struct spart *sparts, size_t nr_parts, size_t nr_gparts,
114
                       size_t nr_sparts, int verbose) {
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

  for (size_t k = 0; k < nr_gparts; ++k) {

    /* We have a DM particle */
    if (gparts[k].type == swift_type_dark_matter) {

      /* Check that it's not linked */
      if (gparts[k].id_or_neg_offset < 0)
        error("DM gpart particle linked to something !");
    }

    /* We have a gas particle */
    else if (gparts[k].type == swift_type_gas) {

      /* Check that it is linked */
      if (gparts[k].id_or_neg_offset > 0)
        error("Gas gpart not linked to anything !");

      /* Find its link */
      const struct part *part = &parts[-gparts[k].id_or_neg_offset];

      /* Check the reverse link */
      if (part->gpart != &gparts[k]) error("Linking problem !");

      /* Check that the particles are at the same place */
      if (gparts[k].x[0] != part->x[0] || gparts[k].x[1] != part->x[1] ||
          gparts[k].x[2] != part->x[2])
142
143
144
145
146
147
        error(
            "Linked particles are not at the same position !\n"
            "gp->x=[%e %e %e] p->x=[%e %e %e] diff=[%e %e %e]",
            gparts[k].x[0], gparts[k].x[1], gparts[k].x[2], part->x[0],
            part->x[1], part->x[2], gparts[k].x[0] - part->x[0],
            gparts[k].x[1] - part->x[1], gparts[k].x[2] - part->x[2]);
148

149
      /* Check that the particles are at the same time */
Matthieu Schaller's avatar
Matthieu Schaller committed
150
      if (gparts[k].time_bin != part->time_bin)
151
        error("Linked particles are not at the same time !");
152
153
154
155
156
157
    }

    else if (gparts[k].type == swift_type_star) {

      /* Check that it is linked */
      if (gparts[k].id_or_neg_offset > 0)
158
        error("Star gpart not linked to anything !");
159
160
161
162
163
164
165
166
167
168

      /* Find its link */
      const struct spart *spart = &sparts[-gparts[k].id_or_neg_offset];

      /* Check the reverse link */
      if (spart->gpart != &gparts[k]) error("Linking problem !");

      /* Check that the particles are at the same place */
      if (gparts[k].x[0] != spart->x[0] || gparts[k].x[1] != spart->x[1] ||
          gparts[k].x[2] != spart->x[2])
169
170
171
172
173
174
        error(
            "Linked particles are not at the same position !\n"
            "gp->x=[%e %e %e] sp->x=[%e %e %e] diff=[%e %e %e]",
            gparts[k].x[0], gparts[k].x[1], gparts[k].x[2], spart->x[0],
            spart->x[1], spart->x[2], gparts[k].x[0] - spart->x[0],
            gparts[k].x[1] - spart->x[1], gparts[k].x[2] - spart->x[2]);
175
176

      /* Check that the particles are at the same time */
Matthieu Schaller's avatar
Matthieu Schaller committed
177
      if (gparts[k].time_bin != spart->time_bin)
178
        error("Linked particles are not at the same time !");
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
    }
  }

  /* Now check that all parts are linked */
  for (size_t k = 0; k < nr_parts; ++k) {

    /* Ok, there is a link */
    if (parts[k].gpart != NULL) {

      /* Check the link */
      if (parts[k].gpart->id_or_neg_offset != -(ptrdiff_t)k) {
        error("Linking problem !");
      }

      /* Check that the particles are at the same place */
      if (parts[k].x[0] != parts[k].gpart->x[0] ||
          parts[k].x[1] != parts[k].gpart->x[1] ||
          parts[k].x[2] != parts[k].gpart->x[2])
        error("Linked particles are not at the same position !");
198
199
200

      /* Check that the particles are at the same time */
      if (parts[k].time_bin != parts[k].gpart->time_bin)
Matthieu Schaller's avatar
Matthieu Schaller committed
201
        error("Linked particles are not at the same time !");
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
    }
  }

  /* Now check that all sparts are linked */
  for (size_t k = 0; k < nr_sparts; ++k) {

    /* Ok, there is a link */
    if (sparts[k].gpart != NULL) {

      /* Check the link */
      if (sparts[k].gpart->id_or_neg_offset != -(ptrdiff_t)k) {
        error("Linking problem !");

        /* Check that the particles are at the same place */
        if (sparts[k].x[0] != sparts[k].gpart->x[0] ||
            sparts[k].x[1] != sparts[k].gpart->x[1] ||
            sparts[k].x[2] != sparts[k].gpart->x[2])
          error("Linked particles are not at the same position !");
220

Matthieu Schaller's avatar
Matthieu Schaller committed
221
222
223
        /* Check that the particles are at the same time */
        if (sparts[k].time_bin != sparts[k].gpart->time_bin)
          error("Linked particles are not at the same time !");
224
225
226
      }
    }
  }
227

228
  if (verbose) message("All links OK");
229
230
}

231
#ifdef WITH_MPI
232
233
234
235
/* MPI data type for the particle transfers */
MPI_Datatype part_mpi_type;
MPI_Datatype xpart_mpi_type;
MPI_Datatype gpart_mpi_type;
236
MPI_Datatype spart_mpi_type;
237
MPI_Datatype multipole_mpi_type;
238
239

/**
240
 * @brief Registers MPI particle types.
241
 */
242
void part_create_mpi_types() {
243

244
  /* This is not the recommended way of doing this.
245
246
247
248
     One should define the structure field by field
     But as long as we don't do serialization via MPI-IO
     we don't really care.
     Also we would have to modify this function everytime something
249
     is added to the part structure. */
250
251
252
253
254
  if (MPI_Type_contiguous(sizeof(struct part) / sizeof(unsigned char), MPI_BYTE,
                          &part_mpi_type) != MPI_SUCCESS ||
      MPI_Type_commit(&part_mpi_type) != MPI_SUCCESS) {
    error("Failed to create MPI type for parts.");
  }
255
256
  if (MPI_Type_contiguous(sizeof(struct xpart) / sizeof(unsigned char),
                          MPI_BYTE, &xpart_mpi_type) != MPI_SUCCESS ||
257
258
259
      MPI_Type_commit(&xpart_mpi_type) != MPI_SUCCESS) {
    error("Failed to create MPI type for xparts.");
  }
260
261
  if (MPI_Type_contiguous(sizeof(struct gpart) / sizeof(unsigned char),
                          MPI_BYTE, &gpart_mpi_type) != MPI_SUCCESS ||
262
263
264
      MPI_Type_commit(&gpart_mpi_type) != MPI_SUCCESS) {
    error("Failed to create MPI type for gparts.");
  }
265
266
267
268
269
  if (MPI_Type_contiguous(sizeof(struct spart) / sizeof(unsigned char),
                          MPI_BYTE, &spart_mpi_type) != MPI_SUCCESS ||
      MPI_Type_commit(&spart_mpi_type) != MPI_SUCCESS) {
    error("Failed to create MPI type for sparts.");
  }
270
271
272
273
274
  if (MPI_Type_contiguous(sizeof(struct multipole) / sizeof(unsigned char),
                          MPI_BYTE, &multipole_mpi_type) != MPI_SUCCESS ||
      MPI_Type_commit(&multipole_mpi_type) != MPI_SUCCESS) {
    error("Failed to create MPI type for multipole.");
  }
275
276
}
#endif