part.c 10.3 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 110 111 112 113 114 115 116
/**
 * @brief Re-link both the #part%s and #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 parts The global #part array in which to find the #gpart offsets.
 * @param sparts The global #spart array in which to find the #gpart offsets.
 */
void part_relink_all_parts_to_gparts(struct gpart *gparts, size_t N,
                                     struct part *parts, struct spart *sparts) {
  for (size_t k = 0; k < N; k++) {
    if (gparts[k].type == swift_type_gas) {
      parts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
    } else if (gparts[k].type == swift_type_star) {
      sparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
    }
  }
}

117 118 119 120 121 122 123 124 125 126 127 128 129
/**
 * @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.
130
 * @param verbose Do we report verbosely in case of success ?
131 132 133
 */
void part_verify_links(struct part *parts, struct gpart *gparts,
                       struct spart *sparts, size_t nr_parts, size_t nr_gparts,
134
                       size_t nr_sparts, int verbose) {
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 */
142
      if (gparts[k].id_or_neg_offset <= 0)
143 144 145 146 147 148 149 150
        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)
151
        error("Gas gpart not linked to anything!");
152 153 154 155 156

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

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

      /* 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])
162
        error(
163
            "Linked particles are not at the same position!\n"
164 165 166 167
            "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]);
168

169
      /* Check that the particles are at the same time */
Matthieu Schaller's avatar
Matthieu Schaller committed
170
      if (gparts[k].time_bin != part->time_bin)
171
        error("Linked particles are not at the same time !");
172 173 174 175 176 177
    }

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

      /* Check that it is linked */
      if (gparts[k].id_or_neg_offset > 0)
178
        error("Star gpart not linked to anything !");
179 180 181 182 183 184 185 186 187 188

      /* 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])
189 190 191 192 193 194
        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]);
195 196

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

  /* 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 !");
218 219 220

      /* 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
221
        error("Linked particles are not at the same time !");
222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239
    }
  }

  /* 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 !");
240

Matthieu Schaller's avatar
Matthieu Schaller committed
241 242 243
        /* 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 !");
244 245 246
      }
    }
  }
247

248
  if (verbose) message("All links OK");
249 250
}

251
#ifdef WITH_MPI
252 253 254 255
/* MPI data type for the particle transfers */
MPI_Datatype part_mpi_type;
MPI_Datatype xpart_mpi_type;
MPI_Datatype gpart_mpi_type;
256
MPI_Datatype spart_mpi_type;
257
MPI_Datatype multipole_mpi_type;
258 259

/**
260
 * @brief Registers MPI particle types.
261
 */
262
void part_create_mpi_types() {
263

264
  /* This is not the recommended way of doing this.
265 266 267 268
     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
269
     is added to the part structure. */
270 271 272 273 274
  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.");
  }
275 276
  if (MPI_Type_contiguous(sizeof(struct xpart) / sizeof(unsigned char),
                          MPI_BYTE, &xpart_mpi_type) != MPI_SUCCESS ||
277 278 279
      MPI_Type_commit(&xpart_mpi_type) != MPI_SUCCESS) {
    error("Failed to create MPI type for xparts.");
  }
280 281
  if (MPI_Type_contiguous(sizeof(struct gpart) / sizeof(unsigned char),
                          MPI_BYTE, &gpart_mpi_type) != MPI_SUCCESS ||
282 283 284
      MPI_Type_commit(&gpart_mpi_type) != MPI_SUCCESS) {
    error("Failed to create MPI type for gparts.");
  }
285 286 287 288 289
  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.");
  }
290 291 292
  if (MPI_Type_contiguous(
          sizeof(struct gravity_tensors) / sizeof(unsigned char), MPI_BYTE,
          &multipole_mpi_type) != MPI_SUCCESS ||
293 294 295
      MPI_Type_commit(&multipole_mpi_type) != MPI_SUCCESS) {
    error("Failed to create MPI type for multipole.");
  }
296 297
}
#endif