Skip to content
Snippets Groups Projects
Select Git revision
  • fd3d64a8c86c95d2eec393834f62bc5ed59bf149
  • master default protected
  • darwin/gear_preSN_feedback
  • zoom-missing-rebuild-time
  • reyz/gear_preSN_feedback
  • darwin/sink_tasks_refactored
  • darwin/gear_radiation
  • darwin/gear_preSN_fbk_merge
  • darwin/gear_mechanical_feedback
  • darwin/gear_chemistry_fluxes
  • zoom_merge protected
  • nickishch/MHD_canvas/SwitchesGradientProjection
  • zoom_mpi_redux
  • MHD_canvas protected
  • karapiperis/plasma_beta_rms_in_tensile_instability_correction_taper_function
  • nickishch/MHD_canvas/OWAR_rhob_div_rhob
  • nickishch/MHD_canvas/OWAR_induction_only
  • nickishch/MHD_canvas/OWAR
  • nickishch/MHD_canvas/Switches
  • zoom_truncate_bkg
  • split-space-split
  • v2025.04 protected
  • v2025.01 protected
  • v1.0.0 protected
  • v0.9.0 protected
  • v0.8.5 protected
  • v0.8.4 protected
  • v0.8.3 protected
  • v0.8.2 protected
  • v0.8.1 protected
  • v0.8.0 protected
  • v0.7.0 protected
  • v0.6.0 protected
  • v0.5.0 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0-pre protected
  • v0.1 protected
  • v0.0 protected
40 results

part.c

Blame
  • user avatar
    Matthieu Schaller authored
    a1f2133c
    History
    part.c 9.57 KiB
    /*******************************************************************************
     * 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. */
    #include "error.h"
    #include "multipole.h"
    #include "part.h"
    
    /**
     * @brief Re-link the #gpart%s associated with the list of #part%s.
     *
     * @param parts The list of #part.
     * @param N The number of particles to re-link;
     * @param offset The offset of #part%s relative to the global parts list.
     */
    void part_relink_gparts_to_parts(struct part *parts, size_t N,
                                     ptrdiff_t offset) {
      for (size_t k = 0; k < N; k++) {
        if (parts[k].gpart) {
          parts[k].gpart->id_or_neg_offset = -(k + offset);
        }
      }
    }
    
    /**
     * @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.
     *
     * @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.
     */
    void part_relink_parts_to_gparts(struct gpart *gparts, size_t N,
                                     struct part *parts) {
      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];
        }
      }
    }
    
    /**
     * @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];
        }
      }
    }
    
    /**
     * @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.
     * @param verbose Do we report verbosely in case of success ?
     */
    void part_verify_links(struct part *parts, struct gpart *gparts,
                           struct spart *sparts, size_t nr_parts, size_t nr_gparts,
                           size_t nr_sparts, int verbose) {
    
      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])
            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]);
    
          /* Check that the particles are at the same time */
          if (gparts[k].time_bin != part->time_bin)
            error("Linked particles are not at the same time !");
        }
    
        else if (gparts[k].type == swift_type_star) {
    
          /* Check that it is linked */
          if (gparts[k].id_or_neg_offset > 0)
            error("Star gpart not linked to anything !");
    
          /* 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])
            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]);
    
          /* Check that the particles are at the same time */
          if (gparts[k].time_bin != spart->time_bin)
            error("Linked particles are not at the same time !");
        }
      }
    
      /* 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 !");
    
          /* Check that the particles are at the same time */
          if (parts[k].time_bin != parts[k].gpart->time_bin)
            error("Linked particles are not at the same time !");
        }
      }
    
      /* 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 !");
    
            /* 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 !");
          }
        }
      }
    
      if (verbose) message("All links OK");
    }
    
    #ifdef WITH_MPI
    /* MPI data type for the particle transfers */
    MPI_Datatype part_mpi_type;
    MPI_Datatype xpart_mpi_type;
    MPI_Datatype gpart_mpi_type;
    MPI_Datatype spart_mpi_type;
    MPI_Datatype multipole_mpi_type;
    
    /**
     * @brief Registers MPI particle types.
     */
    void part_create_mpi_types() {
    
      /* This is not the recommended way of doing this.
         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
         is added to the part structure. */
      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.");
      }
      if (MPI_Type_contiguous(sizeof(struct xpart) / sizeof(unsigned char),
                              MPI_BYTE, &xpart_mpi_type) != MPI_SUCCESS ||
          MPI_Type_commit(&xpart_mpi_type) != MPI_SUCCESS) {
        error("Failed to create MPI type for xparts.");
      }
      if (MPI_Type_contiguous(sizeof(struct gpart) / sizeof(unsigned char),
                              MPI_BYTE, &gpart_mpi_type) != MPI_SUCCESS ||
          MPI_Type_commit(&gpart_mpi_type) != MPI_SUCCESS) {
        error("Failed to create MPI type for gparts.");
      }
      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.");
      }
      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.");
      }
    }
    #endif