Merge branch 'multi_softening' into 'master'

Support for multiple softening lengths in the gravity solver

Closes #599

See merge request !884

examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml

@@ -10,8 +10,8 @@ InternalUnitSystem:
 Gravity:
   eta:          0.025               # Constant dimensionless multiplier for time integration.
   theta:        0.7                 # Opening angle (Multipole acceptance criterion).
-  comoving_softening:     0.1       # Comoving softening length (in internal units).
-  max_physical_softening: 0.1       # Physical softening length (in internal units).
+  max_physical_DM_softening:     0.1       # Physical softening length (in internal units).
+  max_physical_baryon_softening: 0.1       # Physical softening length (in internal units).
 
 # Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
 TimeIntegration:

examples/IsolatedGalaxy/IsolatedGalaxy_potential/isolated_galaxy.yml

@@ -8,11 +8,9 @@ InternalUnitSystem:
 
 # Parameters for the self-gravity scheme
 Gravity:
-  mesh_side_length:       32        # Number of cells along each axis for the periodic gravity mesh.
   eta:          0.025               # Constant dimensionless multiplier for time integration.
   theta:        0.7                 # Opening angle (Multipole acceptance criterion).
-  comoving_softening:     0.0300 # Comoving softening length (in internal units).
-  max_physical_softening: 0.0300    # Physical softening length (in internal units).
+  max_physical_baryon_softening: 0.100  # Physical softening length (in internal units).
 
 # Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
 TimeIntegration:
@@ -34,7 +32,7 @@ Statistics:
   time_first:             0.     # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units)
 
 Scheduler:
-  max_top_level_cells:   96
+  max_top_level_cells:   32
 
 # Parameters related to the initial conditions
 InitialConditions:

examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml

@@ -8,11 +8,10 @@ InternalUnitSystem:
 
 # Parameters for the self-gravity scheme
 Gravity:
-  mesh_side_length:       32        # Number of cells along each axis for the periodic gravity mesh.
-  eta:          0.025               # Constant dimensionless multiplier for time integration.
-  theta:        0.7                 # Opening angle (Multipole acceptance criterion).
-  comoving_softening:     0.01      # Comoving softening length (in internal units).
-  max_physical_softening: 0.01      # Physical softening length (in internal units).
+  eta:          0.025                 # Constant dimensionless multiplier for time integration.
+  theta:        0.7                   # Opening angle (Multipole acceptance criterion).
+  max_physical_DM_softening:     0.1  # Physical softening length (in internal units).
+  max_physical_baryon_softening: 0.1  # Physical softening length (in internal units).
 
 # Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
 TimeIntegration:

examples/PMillennium/PMillennium-1536/p-mill-1536.yml

@@ -43,9 +43,9 @@ Statistics:
 Gravity:
   eta:                    0.025         
   theta:                  0.5          
-  comoving_softening:     0.0208333  # 20.8333 kpc = 1/25 mean inter-particle separation
-  max_physical_softening: 0.0208333  # 20.8333 kpc = 1/25 mean inter-particle separation
-  mesh_side_length:       256
+  comoving_DM_softening:     0.0208333  # 20.8333 kpc = 1/25 mean inter-particle separation
+  max_physical_DM_softening: 0.0208333  # 20.8333 kpc = 1/25 mean inter-particle separation
+  mesh_side_length:       512
   
 # Parameters related to the initial conditions
 InitialConditions:

examples/PMillennium/PMillennium-384/p-mill-384.yml

@@ -43,9 +43,9 @@ Statistics:
 Gravity:
   eta:                    0.025         
   theta:                  0.5          
-  comoving_softening:     0.08333  # 83.333 kpc = 1/25 mean inter-particle separation
-  max_physical_softening: 0.08333  # 83.333 kpc = 1/25 mean inter-particle separation
-  mesh_side_length:       64
+  comoving_DM_softening:     0.08333  # 83.333 kpc = 1/25 mean inter-particle separation
+  max_physical_DM_softening: 0.08333  # 83.333 kpc = 1/25 mean inter-particle separation
+  mesh_side_length:       128
   
 # Parameters related to the initial conditions
 InitialConditions:

examples/PMillennium/PMillennium-768/p-mill-768.yml

@@ -43,9 +43,9 @@ Statistics:
 Gravity:
   eta:                    0.025         
   theta:                  0.5          
-  comoving_softening:     0.041666  # 41.6666 kpc = 1/25 mean inter-particle separation
-  max_physical_softening: 0.041666  # 41.6666 kpc = 1/25 mean inter-particle separation
-  mesh_side_length:       128
+  comoving_DM_softening:     0.041666  # 41.6666 kpc = 1/25 mean inter-particle separation
+  max_physical_DM_softening: 0.041666  # 41.6666 kpc = 1/25 mean inter-particle separation
+  mesh_side_length:       256
   
 # Parameters related to the initial conditions
 InitialConditions:

examples/SantaBarbara/SantaBarbara-128/santa_barbara.yml

@@ -43,9 +43,11 @@ Statistics:
 Gravity:
   eta:                    0.025  
   theta:                  0.5
-  comoving_softening:     0.02    # 20 kpc = 1/25 mean inter-particle separation
-  max_physical_softening: 0.00526 # 20 ckpc = 5.26 pkpc at z=2.8 (EAGLE-like evolution of softening).
-  mesh_side_length:       64
+  comoving_DM_softening:         0.02    # 20 kpc = 1/25 mean inter-particle separation
+  max_physical_DM_softening:     0.00526 # 20 ckpc = 5.26 pkpc at z=2.8 (EAGLE-like evolution of softening).
+  comoving_baryon_softening:     0.02    # 20 kpc = 1/25 mean inter-particle separation
+  max_physical_baryon_softening: 0.00526 # 20 ckpc = 5.26 pkpc at z=2.8 (EAGLE-like evolution of softening).
+  mesh_side_length:       128
 
 # Parameters of the hydro scheme
 SPH:

examples/SantaBarbara/SantaBarbara-256/santa_barbara.yml

@@ -43,9 +43,11 @@ Statistics:
 Gravity:
   eta:                    0.025  
   theta:                  0.5
-  comoving_softening:     0.01    # 10 kpc = 1/25 mean inter-particle separation
-  max_physical_softening: 0.00263 # 10 ckpc = 2.63 pkpc at z=2.8 (EAGLE-like evolution of softening).
-  mesh_side_length:       128
+  comoving_DM_softening:         0.01    # 10 kpc = 1/25 mean inter-particle separation
+  max_physical_DM_softening:     0.00263 # 10 ckpc = 2.63 pkpc at z=2.8 (EAGLE-like evolution of softening).
+  comoving_baryon_softening:     0.01    # 10 kpc = 1/25 mean inter-particle separation
+  max_physical_baryon_softening: 0.00263 # 10 ckpc = 2.63 pkpc at z=2.8 (EAGLE-like evolution of softening).
+  mesh_side_length:       256
 
 # Parameters of the hydro scheme
 SPH:

examples/SmallCosmoVolume/SmallCosmoVolume_DM/small_cosmo_volume_dm.yml

@@ -29,9 +29,9 @@ TimeIntegration:
 # Parameters for the self-gravity scheme
 Gravity:
   eta:          0.025         
-  theta:        0.3           
-  comoving_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
-  max_physical_softening: 0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  theta:        0.5           
+  comoving_DM_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  max_physical_DM_softening: 0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
   mesh_side_length:       64
 
 # Parameters governing the snapshots

examples/SmallCosmoVolume/SmallCosmoVolume_VELOCIraptor/small_cosmo_volume.yml

@@ -22,9 +22,11 @@ TimeIntegration:
 # Parameters for the self-gravity scheme
 Gravity:
   eta:          0.025         
-  theta:        0.3           
-  comoving_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
-  max_physical_softening: 0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  theta:        0.5           
+  comoving_DM_softening:         0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  max_physical_DM_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  comoving_baryon_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  max_physical_baryon_softening: 0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
   mesh_side_length:       64
 
 # Parameters of the hydro scheme

examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml

@@ -23,8 +23,10 @@ TimeIntegration:
 Gravity:
   eta:          0.025         
   theta:        0.3           
-  comoving_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
-  max_physical_softening: 0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  comoving_DM_softening:         0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  max_physical_DM_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  comoving_baryon_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  max_physical_baryon_softening: 0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
   mesh_side_length:       64
 
 # Parameters of the hydro scheme

examples/SmallCosmoVolume/SmallCosmoVolume_hydro/small_cosmo_volume.yml

@@ -22,9 +22,11 @@ TimeIntegration:
 # Parameters for the self-gravity scheme
 Gravity:
   eta:          0.025         
-  theta:        0.3           
-  comoving_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
-  max_physical_softening: 0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  theta:        0.5           
+  comoving_DM_softening:         0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  max_physical_DM_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  comoving_baryon_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  max_physical_baryon_softening: 0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
   mesh_side_length:       64
 
 # Parameters of the hydro scheme

examples/main.c

@@ -440,6 +440,9 @@ int main(int argc, char *argv[]) {
   /* Genesis 1.1: And then, there was time ! */
   clocks_set_cpufreq(cpufreq);
 
+  /* Are we running with gravity */
+  const int with_gravity = (with_self_gravity || with_external_gravity);
+
   /* How vocal are we ? */
   const int talking = (verbose == 1 && myrank == 0) || (verbose == 2);
 
@@ -745,12 +748,6 @@ int main(int argc, char *argv[]) {
     const int generate_gas_in_ics = parser_get_opt_param_int(
         params, "InitialConditions:generate_gas_in_ics", 0);
 
-    /* Some checks that we are not doing something stupid */
-    if (generate_gas_in_ics && flag_entropy_ICs)
-      error("Can't generate gas if the entropy flag is set in the ICs.");
-    if (generate_gas_in_ics && !with_cosmology)
-      error("Can't generate gas if the run is not cosmological.");
-
     /* Initialise the cosmology */
     if (with_cosmology)
       cosmology_init(params, &us, &prog_const, &cosmo);
@@ -811,12 +808,6 @@ int main(int argc, char *argv[]) {
     } else
       bzero(&black_holes_properties, sizeof(struct black_holes_props));
 
-    /* Initialise the gravity properties */
-    bzero(&gravity_properties, sizeof(struct gravity_props));
-    if (with_self_gravity)
-      gravity_props_init(&gravity_properties, params, &prog_const, &cosmo,
-                         with_cosmology, periodic);
-
       /* Initialise the cooling function properties */
 #ifdef COOLING_NONE
     if (with_cooling || with_temperature) {
@@ -868,32 +859,33 @@ int main(int argc, char *argv[]) {
     fflush(stdout);
 
     /* Get ready to read particles of all kinds */
-    size_t Ngas = 0, Ngpart = 0, Nspart = 0, Nbpart = 0;
+    size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0;
     double dim[3] = {0., 0., 0.};
+
     if (myrank == 0) clocks_gettime(&tic);
 #if defined(HAVE_HDF5)
 #if defined(WITH_MPI)
 #if defined(HAVE_PARALLEL_HDF5)
     read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
-                     &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
-                     with_hydro, (with_external_gravity || with_self_gravity),
-                     with_stars, with_black_holes, cleanup_h, cleanup_sqrt_a,
-                     cosmo.h, cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD,
-                     MPI_INFO_NULL, nr_threads, dry_run);
+                     &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
+                     &flag_entropy_ICs, with_hydro, with_gravity, with_stars,
+                     with_black_holes, cleanup_h, cleanup_sqrt_a, cosmo.h,
+                     cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL,
+                     nr_threads, dry_run);
 #else
     read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
-                   &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
-                   with_hydro, (with_external_gravity || with_self_gravity),
-                   with_stars, with_black_holes, cleanup_h, cleanup_sqrt_a,
-                   cosmo.h, cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD,
-                   MPI_INFO_NULL, nr_threads, dry_run);
+                   &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
+                   &flag_entropy_ICs, with_hydro, with_gravity, with_stars,
+                   with_black_holes, cleanup_h, cleanup_sqrt_a, cosmo.h,
+                   cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL,
+                   nr_threads, dry_run);
 #endif
 #else
     read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
-                   &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
-                   with_hydro, (with_external_gravity || with_self_gravity),
-                   with_stars, with_black_holes, cleanup_h, cleanup_sqrt_a,
-                   cosmo.h, cosmo.a, nr_threads, dry_run);
+                   &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
+                   &flag_entropy_ICs, with_hydro, with_gravity, with_stars,
+                   with_black_holes, cleanup_h, cleanup_sqrt_a, cosmo.h,
+                   cosmo.a, nr_threads, dry_run);
 #endif
 #endif
     if (myrank == 0) {
@@ -903,6 +895,12 @@ int main(int argc, char *argv[]) {
       fflush(stdout);
     }
 
+    /* Some checks that we are not doing something stupid */
+    if (generate_gas_in_ics && flag_entropy_ICs)
+      error("Can't generate gas if the entropy flag is set in the ICs.");
+    if (generate_gas_in_ics && !with_cosmology)
+      error("Can't generate gas if the run is not cosmological.");
+
 #ifdef SWIFT_DEBUG_CHECKS
     /* Check once and for all that we don't have unwanted links */
     if (!with_stars && !dry_run) {
@@ -925,23 +923,46 @@ int main(int argc, char *argv[]) {
 #endif
 
     /* Get the total number of particles across all nodes. */
-    long long N_total[4] = {0, 0, 0, 0};
+    long long N_total[swift_type_count + 1] = {0};
+    long long Nbaryons = Ngas + Nspart + Nbpart;
 #if defined(WITH_MPI)
-    long long N_long[4] = {Ngas, Ngpart, Nspart, Nbpart};
-    MPI_Allreduce(&N_long, &N_total, 4, MPI_LONG_LONG_INT, MPI_SUM,
-                  MPI_COMM_WORLD);
+    long long N_long[swift_type_count + 1] = {0};
+    N_long[swift_type_gas] = Ngas;
+    N_long[swift_type_dark_matter] =
+        with_gravity ? Ngpart - Ngpart_background - Nbaryons : 0;
+    N_long[swift_type_dark_matter_background] = Ngpart_background;
+    N_long[swift_type_stars] = Nspart;
+    N_long[swift_type_black_hole] = Nbpart;
+    N_long[swift_type_count] = Ngpart;
+    MPI_Allreduce(&N_long, &N_total, swift_type_count + 1, MPI_LONG_LONG_INT,
+                  MPI_SUM, MPI_COMM_WORLD);
 #else
-    N_total[0] = Ngas;
-    N_total[1] = Ngpart;
-    N_total[2] = Nspart;
-    N_total[3] = Nbpart;
+    N_total[swift_type_gas] = Ngas;
+    N_total[swift_type_dark_matter] =
+        with_gravity ? Ngpart - Ngpart_background - Nbaryons : 0;
+    N_total[swift_type_dark_matter_background] = Ngpart_background;
+    N_total[swift_type_stars] = Nspart;
+    N_total[swift_type_black_hole] = Nbpart;
+    N_total[swift_type_count] = Ngpart;
 #endif
 
     if (myrank == 0)
       message(
           "Read %lld gas particles, %lld stars particles, %lld black hole "
-          "particles and %lld gparts from the ICs.",
-          N_total[0], N_total[2], N_total[3], N_total[1]);
+          "particles, %lld DM particles and %lld DM background particles from "
+          "the ICs.",
+          N_total[swift_type_gas], N_total[swift_type_stars],
+          N_total[swift_type_black_hole], N_total[swift_type_dark_matter],
+          N_total[swift_type_dark_matter_background]);
+
+    const int with_DM_particles = N_total[swift_type_dark_matter] > 0;
+    const int with_baryon_particles =
+        (N_total[swift_type_gas] + N_total[swift_type_stars] +
+         N_total[swift_type_black_hole]) > 0;
+
+    /* Do we have background DM particles? */
+    const int with_DM_background_particles =
+        N_total[swift_type_dark_matter_background] > 0;
 
     /* Verify that the fields to dump actually exist */
     if (myrank == 0) io_check_output_fields(params, N_total);
@@ -950,8 +971,8 @@ int main(int argc, char *argv[]) {
     if (myrank == 0) clocks_gettime(&tic);
     space_init(&s, params, &cosmo, dim, parts, gparts, sparts, bparts, Ngas,
                Ngpart, Nspart, Nbpart, periodic, replicate, generate_gas_in_ics,
-               with_hydro, with_self_gravity, with_star_formation, talking,
-               dry_run);
+               with_hydro, with_self_gravity, with_star_formation,
+               with_DM_background_particles, talking, dry_run);
 
     if (myrank == 0) {
       clocks_gettime(&toc);
@@ -960,6 +981,14 @@ int main(int argc, char *argv[]) {
       fflush(stdout);
     }
 
+    /* Initialise the gravity properties */
+    bzero(&gravity_properties, sizeof(struct gravity_props));
+    if (with_self_gravity)
+      gravity_props_init(&gravity_properties, params, &prog_const, &cosmo,
+                         with_cosmology, with_baryon_particles,
+                         with_DM_particles, with_DM_background_particles,
+                         periodic);
+
     /* Initialise the external potential properties */
     bzero(&potential, sizeof(struct external_potential));
     if (with_external_gravity)
@@ -984,19 +1013,24 @@ int main(int argc, char *argv[]) {
     if (with_cosmology && with_self_gravity && !dry_run)
       space_check_cosmology(&s, &cosmo, myrank);
 
-/* Also update the total counts (in case of changes due to replication) */
+    /* Also update the total counts (in case of changes due to replication) */
+    Nbaryons = s.nr_parts + s.nr_sparts + s.nr_bparts;
 #if defined(WITH_MPI)
-    N_long[0] = s.nr_parts;
-    N_long[1] = s.nr_gparts;
-    N_long[2] = s.nr_sparts;
-    N_long[3] = s.nr_bparts;
-    MPI_Allreduce(&N_long, &N_total, 4, MPI_LONG_LONG_INT, MPI_SUM,
-                  MPI_COMM_WORLD);
+    N_long[swift_type_gas] = s.nr_parts;
+    N_long[swift_type_dark_matter] =
+        with_gravity ? s.nr_gparts - Ngpart_background - Nbaryons : 0;
+    N_long[swift_type_count] = s.nr_gparts;
+    N_long[swift_type_stars] = s.nr_sparts;
+    N_long[swift_type_black_hole] = s.nr_bparts;
+    MPI_Allreduce(&N_long, &N_total, swift_type_count + 1, MPI_LONG_LONG_INT,
+                  MPI_SUM, MPI_COMM_WORLD);
 #else
-    N_total[0] = s.nr_parts;
-    N_total[1] = s.nr_gparts;
-    N_total[2] = s.nr_sparts;
-    N_total[3] = s.nr_bparts;
+    N_total[swift_type_gas] = s.nr_parts;
+    N_total[swift_type_dark_matter] =
+        with_gravity ? s.nr_gparts - Ngpart_background - Nbaryons : 0;
+    N_total[swift_type_count] = s.nr_gparts;
+    N_total[swift_type_stars] = s.nr_sparts;
+    N_total[swift_type_black_hole] = s.nr_bparts;
 #endif
 
     /* Say a few nice things about the space we just created. */
@@ -1020,7 +1054,7 @@ int main(int argc, char *argv[]) {
           "ERROR: Running with hydrodynamics but no gas particles found in the "
           "ICs!");
     }
-    if ((with_self_gravity || with_external_gravity) && N_total[1] == 0) {
+    if (with_gravity && N_total[1] == 0) {
       error(
           "ERROR: Running with gravity but no gravity particles found in "
           "the ICs!");
@@ -1063,12 +1097,14 @@ int main(int argc, char *argv[]) {
 
     /* Initialize the engine with the space and policies. */
     if (myrank == 0) clocks_gettime(&tic);
-    engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2], N_total[3],
-                engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
-                &hydro_properties, &entropy_floor, &gravity_properties,
-                &stars_properties, &black_holes_properties,
-                &feedback_properties, &mesh, &potential, &cooling_func,
-                &starform, &chemistry, &fof_properties);
+    engine_init(
+        &e, &s, params, N_total[swift_type_gas], N_total[swift_type_count],
+        N_total[swift_type_stars], N_total[swift_type_black_hole],
+        N_total[swift_type_dark_matter_background], engine_policies, talking,
+        &reparttype, &us, &prog_const, &cosmo, &hydro_properties,
+        &entropy_floor, &gravity_properties, &stars_properties,
+        &black_holes_properties, &feedback_properties, &mesh, &potential,
+        &cooling_func, &starform, &chemistry, &fof_properties);
     engine_config(/*restart=*/0, /*fof=*/0, &e, params, nr_nodes, myrank,
                   nr_threads, with_aff, talking, restart_file);
 
@@ -1081,12 +1117,13 @@ int main(int argc, char *argv[]) {
 
     /* Get some info to the user. */
     if (myrank == 0) {
-      long long N_DM = N_total[1] - N_total[2] - N_total[3] - N_total[0];
+      const long long N_DM = N_total[swift_type_dark_matter] +
+                             N_total[swift_type_dark_matter_background];
       message(
           "Running on %lld gas particles, %lld stars particles %lld black "
           "hole particles and %lld DM particles (%lld gravity particles)",
-          N_total[0], N_total[2], N_total[3], N_total[1] > 0 ? N_DM : 0,
-          N_total[1]);
+          N_total[swift_type_gas], N_total[swift_type_stars],
+          N_total[swift_type_black_hole], N_DM, N_total[swift_type_count]);
       message(
           "from t=%.3e until t=%.3e with %d ranks, %d threads / rank and %d "
           "task queues / rank (dt_min=%.3e, dt_max=%.3e)...",

examples/main_fof.c

@@ -409,11 +409,6 @@ int main(int argc, char *argv[]) {
   /* Initialise the cosmology */
   cosmology_init(params, &us, &prog_const, &cosmo);
 
-  /* Initialise the gravity scheme */
-  bzero(&gravity_properties, sizeof(struct gravity_props));
-  gravity_props_init(&gravity_properties, params, &prog_const, &cosmo,
-                     /*with_cosmology=*/1, periodic);
-
   /* Initialise the FOF properties */
   bzero(&fof_properties, sizeof(struct fof_props));
   if (with_fof) fof_init(&fof_properties, params, &prog_const, &us);
@@ -428,28 +423,32 @@ int main(int argc, char *argv[]) {
 
   /* Get ready to read particles of all kinds */
   int flag_entropy_ICs = 0;
-  size_t Ngas = 0, Ngpart = 0, Nspart = 0, Nbpart = 0;
+  size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0;
   double dim[3] = {0., 0., 0.};
+
   if (myrank == 0) clocks_gettime(&tic);
 #if defined(HAVE_HDF5)
 #if defined(WITH_MPI)
 #if defined(HAVE_PARALLEL_HDF5)
   read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
-                   &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
-                   with_hydro, /*with_grav=*/1, with_stars, with_black_holes,
-                   cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
-                   nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
+                   &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
+                   &flag_entropy_ICs, with_hydro,
+                   /*with_grav=*/1, with_stars, with_black_holes, cleanup_h,
+                   cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes,
+                   MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
                    /*dry_run=*/0);
 #else
   read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas,
-                 &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro,
+                 &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
+                 &flag_entropy_ICs, with_hydro,
                  /*with_grav=*/1, with_stars, with_black_holes, cleanup_h,
                  cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes,
                  MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, /*dry_run=*/0);
 #endif
 #else
   read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas,
-                 &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro,
+                 &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
+                 &flag_entropy_ICs, with_hydro,
                  /*with_grav=*/1, with_stars, with_black_holes, cleanup_h,
                  cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads, /*dry_run=*/0);
 #endif
@@ -468,30 +467,51 @@ int main(int argc, char *argv[]) {
 #endif
 
   /* Get the total number of particles across all nodes. */
-  long long N_total[4] = {0, 0, 0};
+  long long N_total[swift_type_count + 1] = {0};
+  long long Nbaryons = Ngas + Nspart + Nbpart;
 #if defined(WITH_MPI)
-  long long N_long[4] = {Ngas, Ngpart, Nspart, Nbpart};
-  MPI_Allreduce(&N_long, &N_total, 4, MPI_LONG_LONG_INT, MPI_SUM,
-                MPI_COMM_WORLD);
+  long long N_long[swift_type_count + 1] = {0};
+  N_long[swift_type_gas] = Ngas;
+  N_long[swift_type_dark_matter] = Ngpart - Ngpart_background - Nbaryons;
+  N_long[swift_type_dark_matter_background] = Ngpart_background;
+  N_long[swift_type_stars] = Nspart;
+  N_long[swift_type_black_hole] = Nbpart;
+  N_long[swift_type_count] = Ngpart;
+  MPI_Allreduce(&N_long, &N_total, swift_type_count + 1, MPI_LONG_LONG_INT,
+                MPI_SUM, MPI_COMM_WORLD);
 #else
-  N_total[0] = Ngas;
-  N_total[1] = Ngpart;
-  N_total[2] = Nspart;
-  N_total[3] = Nbpart;
+  N_total[swift_type_gas] = Ngas;
+  N_total[swift_type_dark_matter] = Ngpart - Ngpart_background - Nbaryons;
+  N_total[swift_type_dark_matter_background] = Ngpart_background;
+  N_total[swift_type_stars] = Nspart;
+  N_total[swift_type_black_hole] = Nbpart;
+  N_total[swift_type_count] = Ngpart;
 #endif
 
   if (myrank == 0)
     message(
         "Read %lld gas particles, %lld stars particles, %lld black hole "
-        "particles and %lld gparts from the ICs.",
-        N_total[0], N_total[2], N_total[3], N_total[1]);
+        "particles, %lld DM particles and %lld DM background particles from "
+        "the ICs.",
+        N_total[swift_type_gas], N_total[swift_type_stars],
+        N_total[swift_type_black_hole], N_total[swift_type_dark_matter],
+        N_total[swift_type_dark_matter_background]);
+
+  const int with_DM_particles = N_total[swift_type_dark_matter] > 0;
+  const int with_baryon_particles =
+      (N_total[swift_type_gas] + N_total[swift_type_stars] +
+       N_total[swift_type_black_hole]) > 0;
+
+  /* Do we have background DM particles? */
+  const int with_DM_background_particles =
+      N_total[swift_type_dark_matter_background] > 0;
 
   /* Initialize the space with these data. */
   if (myrank == 0) clocks_gettime(&tic);
   space_init(&s, params, &cosmo, dim, parts, gparts, sparts, bparts, Ngas,
              Ngpart, Nspart, Nbpart, periodic, replicate,
              /*generate_gas_in_ics=*/0, /*hydro=*/N_total[0] > 0, /*gravity=*/1,
-             /*with_star_formation=*/0, talking,
+             /*with_star_formation=*/0, with_DM_background_particles, talking,
              /*dry_run=*/0);
 
   if (myrank == 0) {
@@ -501,6 +521,12 @@ int main(int argc, char *argv[]) {
     fflush(stdout);
   }
 
+  /* Initialise the gravity scheme */
+  bzero(&gravity_properties, sizeof(struct gravity_props));
+  gravity_props_init(&gravity_properties, params, &prog_const, &cosmo,
+                     /*with_cosmology=*/1, with_baryon_particles,
+                     with_DM_particles, with_DM_background_particles, periodic);
+
   /* Initialise the long-range gravity mesh */
   if (periodic) {
 #ifdef HAVE_FFTW
@@ -513,17 +539,22 @@ int main(int argc, char *argv[]) {
     pm_mesh_init_no_mesh(&mesh, s.dim);
   }
 
-    /* Also update the total counts (in case of changes due to replication) */
+  /* Also update the total counts (in case of changes due to replication) */
+  Nbaryons = s.nr_parts + s.nr_sparts + s.nr_bparts;
 #if defined(WITH_MPI)
-  N_long[0] = s.nr_parts;
-  N_long[1] = s.nr_gparts;
-  N_long[2] = s.nr_sparts;
-  MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
-                MPI_COMM_WORLD);
+  N_long[swift_type_gas] = s.nr_parts;
+  N_long[swift_type_dark_matter] = s.nr_gparts - Ngpart_background - Nbaryons;
+  N_long[swift_type_count] = s.nr_gparts;
+  N_long[swift_type_stars] = s.nr_sparts;
+  N_long[swift_type_black_hole] = s.nr_bparts;
+  MPI_Allreduce(&N_long, &N_total, swift_type_count + 1, MPI_LONG_LONG_INT,
+                MPI_SUM, MPI_COMM_WORLD);
 #else
-  N_total[0] = s.nr_parts;
-  N_total[1] = s.nr_gparts;
-  N_total[2] = s.nr_sparts;
+  N_total[swift_type_gas] = s.nr_parts;
+  N_total[swift_type_dark_matter] = s.nr_gparts - Ngpart_background - Nbaryons;
+  N_total[swift_type_count] = s.nr_gparts;
+  N_total[swift_type_stars] = s.nr_sparts;
+  N_total[swift_type_black_hole] = s.nr_bparts;
 #endif
 
   /* Say a few nice things about the space we just created. */
@@ -548,14 +579,16 @@ int main(int argc, char *argv[]) {
 
   /* Initialize the engine with the space and policies. */
   if (myrank == 0) clocks_gettime(&tic);
-  engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2], N_total[3],
-              engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
-              /*hydro_properties=*/NULL, /*entropy_floor=*/NULL,
-              &gravity_properties,
-              /*stars_properties=*/NULL, /*black_holes_properties=*/NULL,
-              /*feedback_properties=*/NULL, &mesh, /*potential=*/NULL,
-              /*cooling_func=*/NULL,
-              /*starform=*/NULL, /*chemistry=*/NULL, &fof_properties);
+  engine_init(
+      &e, &s, params, N_total[swift_type_gas], N_total[swift_type_count],
+      N_total[swift_type_stars], N_total[swift_type_black_hole],
+      N_total[swift_type_dark_matter_background], engine_policies, talking,
+      &reparttype, &us, &prog_const, &cosmo,
+      /*hydro_properties=*/NULL, /*entropy_floor=*/NULL, &gravity_properties,
+      /*stars_properties=*/NULL, /*black_holes_properties=*/NULL,
+      /*feedback_properties=*/NULL, &mesh, /*potential=*/NULL,
+      /*cooling_func=*/NULL,
+      /*starform=*/NULL, /*chemistry=*/NULL, &fof_properties);
   engine_config(/*restart=*/0, /*fof=*/1, &e, params, nr_nodes, myrank,
                 nr_threads, with_aff, talking, NULL);
 
@@ -568,12 +601,13 @@ int main(int argc, char *argv[]) {
 
   /* Get some info to the user. */
   if (myrank == 0) {
-    long long N_DM = N_total[1] - N_total[2] - N_total[3] - N_total[0];
+    const long long N_DM = N_total[swift_type_dark_matter] +
+                           N_total[swift_type_dark_matter_background];
     message(
-        "Running on %lld gas particles, %lld stars particles %lld black "
+        "Running FOF on %lld gas particles, %lld stars particles %lld black "
         "hole particles and %lld DM particles (%lld gravity particles)",
-        N_total[0], N_total[2], N_total[3], N_total[1] > 0 ? N_DM : 0,
-        N_total[1]);
+        N_total[swift_type_gas], N_total[swift_type_stars],
+        N_total[swift_type_black_hole], N_DM, N_total[swift_type_count]);
     message(
         "from t=%.3e until t=%.3e with %d ranks, %d threads / rank and %d "
         "task queues / rank (dt_min=%.3e, dt_max=%.3e)...",

examples/nIFTyCluster/Baryonic/nifty.yml

@@ -37,13 +37,14 @@ Statistics:
 
 # Parameters for the self-gravity scheme
 Gravity:
-  eta:                    0.025    # Constant dimensionless multiplier for time integration. 
-  theta:                  0.7     # Opening angle (Multipole acceptance criterion)
-  comoving_softening:     20.0 # Comoving softening length (in internal units).
-  max_physical_softening: 5.0    # Physical softening length (in internal units).
-  mesh_side_length:       512
-  softening_ratio:        0.04   # 1/25
-  softening_ratio_background:        0.04   # 1/25
+  eta:                           0.025  # Constant dimensionless multiplier for time integration. 
+  theta:                         0.5    # Opening angle (Multipole acceptance criterion)
+  comoving_DM_softening:        20.0    # Comoving softening length (in internal units).
+  comoving_baryon_softening:    20.0    # Comoving softening length (in internal units).
+  max_physical_DM_softening:     5.0    # Max physical softening length (in internal units).
+  max_physical_baryon_softening: 5.0    # Max physical softening length (in internal units).
+  softening_ratio_background:    0.04   # 1/25th of mean inter-particle separation
+  mesh_side_length:              512
 
 # Parameters for the hydrodynamics scheme
 SPH: