Merge branch 'parallel_mesh_assignment' into 'master'

Parallel mesh assignment also in single-node case

See merge request !1509

doc/RTD/source/ParameterFiles/parameter_description.rst

@@ -293,6 +293,9 @@ Particle-Mesh part of the calculation. The last five are optional:
 
 * The number cells along each axis of the mesh :math:`N`: ``mesh_side_length``,
 * Whether or not to use a distributed mesh when running over MPI: ``distributed_mesh`` (default: ``0``),
+* Whether or not to use local patches instead of direct atomic operations to
+  write to the mesh in the non-MPI case (this is a performance tuning
+  parameter): ``mesh_uses_local_patches`` (default: ``1``),
 * The mesh smoothing scale in units of the mesh cell-size :math:`a_{\rm
   smooth}`: ``a_smooth`` (default: ``1.25``),
 * The scale above which the short-range forces are assumed to be 0 (in units of

examples/parameter_example.yml

@@ -79,6 +79,7 @@ Stars:
 Gravity:
   mesh_side_length:              128       # Number of cells along each axis for the periodic gravity mesh (must be even).
   distributed_mesh:              0         # (Optional) Are we using a distributed mesh when running over MPI (necessary for meshes > 1290^3)
+  mesh_uses_local_patches:       1         # (Optional) Are we using thread-local patches (1) or direct atomic writes to the global mesh (0) in the non-MPI case?
   eta:                           0.025     # Constant dimensionless multiplier for time integration.
   MAC:                           adaptive  # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'.
   epsilon_fmm:                   0.001     # Tolerance parameter for the adaptive multipole acceptance criterion.

src/gravity_properties.c

@@ -64,6 +64,8 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
     p->distributed_mesh =
         parser_get_opt_param_int(params, "Gravity:distributed_mesh",
                                  gravity_props_default_distributed_mesh);
+    p->mesh_uses_local_patches =
+        parser_get_opt_param_int(params, "Gravity:mesh_uses_local_patches", 1);
     p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth",
                                              gravity_props_default_a_smooth);
     p->r_cut_max_ratio = parser_get_opt_param_float(

src/gravity_properties.h

@@ -121,6 +121,10 @@ struct gravity_props {
   /*! Whether mesh is distributed between MPI ranks when we use MPI  */
   int distributed_mesh;
 
+  /*! Whether or not to use local patches rather than
+   * direct atomic writes to the mesh when running without MPI */
+  int mesh_uses_local_patches;
+
   /*! Mesh smoothing scale in units of top-level cell size */
   float a_smooth;
 

src/mesh_gravity.c

@@ -215,6 +215,7 @@ struct cic_mapper_data {
   double* rho;
   double* potential;
   int N;
+  int use_local_patches;
   double fac;
   double dim[3];
   float const_G;
@@ -260,10 +261,8 @@ void cell_gpart_to_mesh_CIC_mapper(void* map_data, int num, void* extra) {
   /* Pointer to the chunk to be processed */
   int* local_cells = (int*)map_data;
 
-  // MATTHIEU: This could in principle be improved by creating a local mesh
-  //           with just the extent required for the cell. Assignment can
-  //           then be done without atomics. That local mesh is then added
-  //           atomically to the global one.
+  /* A temporary patch of the global mesh */
+  struct pm_mesh_patch patch;
 
   /* Loop over the elements assigned to this thread */
   for (int i = 0; i < num; ++i) {
@@ -271,8 +270,26 @@ void cell_gpart_to_mesh_CIC_mapper(void* map_data, int num, void* extra) {
     /* Pointer to local cell */
     const struct cell* c = &cells[local_cells[i]];
 
-    /* Assign this cell's content to the mesh */
-    cell_gpart_to_mesh_CIC(c, rho, N, fac, dim, nu_model);
+    /* Skip empty cells */
+    if (c->grav.count == 0) continue;
+
+    if (data->use_local_patches) {
+
+      /* Do a CIC interpolation of all the particles in this cell onto
+         the local patch (allocates memory in the patch) */
+      accumulate_cell_to_local_patch(N, fac, dim, c, &patch, nu_model);
+
+      /* Copy the local patch values back onto the global mesh */
+      pm_add_patch_to_global_mesh(rho, &patch);
+
+      /* Free the allocated memory */
+      pm_mesh_patch_clean(&patch);
+
+    } else {
+
+      /* Assign this cell's content directly atomically to the mesh */
+      cell_gpart_to_mesh_CIC(c, rho, N, fac, dim, nu_model);
+    }
   }
 }
 
@@ -880,6 +897,7 @@ void compute_potential_global(struct pm_mesh* mesh, const struct space* s,
   data.rho = rho;
   data.potential = NULL;
   data.N = N;
+  data.use_local_patches = mesh->use_local_patches;
   data.fac = cell_fac;
   data.dim[0] = dim[0];
   data.dim[1] = dim[1];
@@ -1136,6 +1154,7 @@ void pm_mesh_init(struct pm_mesh* mesh, const struct gravity_props* props,
   mesh->periodic = 1;
   mesh->N = N;
   mesh->distributed_mesh = props->distributed_mesh;
+  mesh->use_local_patches = props->mesh_uses_local_patches;
   mesh->dim[0] = dim[0];
   mesh->dim[1] = dim[1];
   mesh->dim[2] = dim[2];

src/mesh_gravity.h

@@ -50,6 +50,10 @@ struct pm_mesh {
   /*! Whether mesh is distributed between MPI ranks */
   int distributed_mesh;
 
+  /*! Whether or not to use local patches rather than
+   * direct atomic writes to the mesh when running without MPI */
+  int use_local_patches;
+
   /*! Integer time-step end of the mesh force for the last step */
   integertime_t ti_end_mesh_last;
 

src/mesh_gravity_mpi.c

@@ -161,6 +161,9 @@ void accumulate_cell_to_local_patches_mapper(void *map_data, int num,
     /* Pointer to local cell */
     const struct cell *c = &cells[local_cells[i]];
 
+    /* Skip empty cells */
+    if (c->grav.count == 0) continue;
+
     /* Assign this cell's content to the mesh */
     accumulate_cell_to_local_patch(N, fac, dim, c, &local_patches[i], nu_model);
   }

src/mesh_gravity_mpi.h

@@ -28,6 +28,12 @@ struct cell;
 struct threadpool;
 struct pm_mesh;
 struct pm_mesh_patch;
+struct neutrino_model;
+
+void accumulate_cell_to_local_patch(const int N, const double fac,
+                                    const double *dim, const struct cell *cell,
+                                    struct pm_mesh_patch *patch,
+                                    const struct neutrino_model *nu_model);
 
 void mpi_mesh_accumulate_gparts_to_local_patches(
     struct threadpool *tp, const int N, const double fac, const struct space *s,

src/mesh_gravity_patch.c

@@ -26,6 +26,7 @@
 #include "mesh_gravity_patch.h"
 
 /* Local includes. */
+#include "atomic.h"
 #include "cell.h"
 #include "error.h"
 #include "row_major_id.h"
@@ -99,6 +100,45 @@ void pm_mesh_patch_init(struct pm_mesh_patch *patch, const struct cell *cell,
     error("Failed to allocate array for mesh patch!");
 }
 
+/**
+ * @brief Write the content of a mesh patch back to the global mesh
+ * using atomic operations.
+ *
+ * @param global_mesh The global mesh to write to.
+ * @param patch The #pm_mesh_patch object to write from.
+ */
+void pm_add_patch_to_global_mesh(double *const global_mesh,
+                                 const struct pm_mesh_patch *patch) {
+
+  const int N = patch->N;
+  const int size_i = patch->mesh_size[0];
+  const int size_j = patch->mesh_size[1];
+  const int size_k = patch->mesh_size[2];
+  const int mesh_min_i = patch->mesh_min[0];
+  const int mesh_min_j = patch->mesh_min[1];
+  const int mesh_min_k = patch->mesh_min[2];
+
+  /* Remind the compiler that the arrays are nicely aligned */
+  swift_declare_aligned_ptr(const double, mesh, patch->mesh,
+                            SWIFT_CACHE_ALIGNMENT);
+
+  for (int i = 0; i < size_i; ++i) {
+    for (int j = 0; j < size_j; ++j) {
+      for (int k = 0; k < size_k; ++k) {
+
+        const int ii = i + mesh_min_i;
+        const int jj = j + mesh_min_j;
+        const int kk = k + mesh_min_k;
+
+        const int patch_index = pm_mesh_patch_index(patch, i, j, k);
+        const int mesh_index = row_major_id_periodic(ii, jj, kk, N);
+
+        atomic_add_d(&global_mesh[mesh_index], mesh[patch_index]);
+      }
+    }
+  }
+}
+
 /**
  * @brief Set all values in a mesh patch to zero
  *
@@ -106,9 +146,12 @@ void pm_mesh_patch_init(struct pm_mesh_patch *patch, const struct cell *cell,
  */
 void pm_mesh_patch_zero(struct pm_mesh_patch *patch) {
 
+  /* Remind the compiler that the arrays are nicely aligned */
+  swift_declare_aligned_ptr(double, mesh, patch->mesh, SWIFT_CACHE_ALIGNMENT);
+
   const int num =
       patch->mesh_size[0] * patch->mesh_size[1] * patch->mesh_size[2];
-  memset(patch->mesh, 0, num * sizeof(double));
+  memset(mesh, 0, num * sizeof(double));
 }
 
 /**
@@ -118,7 +161,7 @@ void pm_mesh_patch_zero(struct pm_mesh_patch *patch) {
  */
 void pm_mesh_patch_clean(struct pm_mesh_patch *patch) {
 
-  swift_free("mesh_patch", patch->mesh);
+  if (patch->mesh) swift_free("mesh_patch", patch->mesh);
 
   /* Zero everything and give a silly mesh size to help debugging */
   memset(patch, 0, sizeof(struct pm_mesh_patch));

src/mesh_gravity_patch.h

@@ -167,4 +167,7 @@ __attribute__((always_inline)) INLINE static void pm_mesh_patch_CIC_set(
   mesh[pm_mesh_patch_index(patch, i + 1, j + 1, k + 1)] += value * dx * dy * dz;
 }
 
+void pm_add_patch_to_global_mesh(double *const global_mesh,
+                                 const struct pm_mesh_patch *patch);
+
 #endif