diff --git a/configure.ac b/configure.ac
index 28a74d547225ff010eb29cd43182d57bac699e82..338edec60f956c37f666f7592a931b2c20a9f6e8 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1390,24 +1390,21 @@ AM_CONDITIONAL([HAVESTANDALONEFOF],[test $enable_standalone_fof = "yes"])
# Gravity scheme.
AC_ARG_WITH([gravity],
[AS_HELP_STRING([--with-gravity=<scheme>],
- [Gravity scheme to use @<:@default, with-potential, default: default@:>@]
+ [Gravity scheme to use @<:@basic, with-potential, with-multi-softening default: with-multi-softening@:>@]
)],
[with_gravity="$withval"],
- [with_gravity="default"]
+ [with_gravity="with-multi-softening"]
)
-if test "$with_subgrid" = "EAGLE"; then
- if test "$with_gravity" = "default"; then
- with_gravity="with-potential"
- fi
-fi
-
case "$with_gravity" in
with-potential)
- AC_DEFINE([POTENTIAL_GRAVITY], [1], [Gravity scheme with potential calculation])
+ AC_DEFINE([POTENTIAL_GRAVITY], [1], [Basic gravity scheme with potential calculation])
;;
- default)
- AC_DEFINE([DEFAULT_GRAVITY], [1], [Default gravity scheme])
+ with-multi-softening)
+ AC_DEFINE([MULTI_SOFTENING_GRAVITY], [1], [Gravity scheme with per-particle type softening value and background particles])
+ ;;
+ basic)
+ AC_DEFINE([DEFAULT_GRAVITY], [1], [Basic gravity scheme])
;;
*)
AC_MSG_ERROR([Unknown gravity scheme: $with_gravity])
@@ -1982,13 +1979,6 @@ AC_DEFINE_UNQUOTED([SELF_GRAVITY_MULTIPOLE_ORDER], [$with_multipole_order], [Mul
AC_PATH_PROG([GIT_CMD], [git])
AC_SUBST([GIT_CMD])
-# Check that the gravity model is compatible with the subgrid
-if test $with_black_holes = "EAGLE"; then
- if test $with_gravity != "with-potential"; then
- AC_MSG_ERROR([The EAGLE BH model needs the gravity scheme to provide potentials. The code must be compile with --with-gravity=with-potential.])
- fi
-fi
-
# Make the documentation. Add conditional to handle disable option.
DX_INIT_DOXYGEN(libswift,doc/Doxyfile,doc/)
AM_CONDITIONAL([HAVE_DOXYGEN], [test "$ac_cv_path_ac_pt_DX_DOXYGEN" != ""])
diff --git a/doc/RTD/source/InitialConditions/index.rst b/doc/RTD/source/InitialConditions/index.rst
index e585c9aa55f269ebbbf9b2d83034b96a688a99f4..4532fdde8818088cf89cee5fac7cacfbee2dd67b 100644
--- a/doc/RTD/source/InitialConditions/index.rst
+++ b/doc/RTD/source/InitialConditions/index.rst
@@ -44,27 +44,29 @@ There are several groups that contain 'auxiliary' information, such as
the particles. Some types are currently ignored by SWIFT but are kept in the
file format for compatibility reasons.
-+---------------------+------------------------+----------------------------+
-| HDF5 Group Name | Physical Particle Type | In code ``enum part_type`` |
-+=====================+========================+============================+
-| ``/PartType0/`` | Gas | ``swift_type_gas`` |
-+---------------------+------------------------+----------------------------+
-| ``/PartType1/`` | Dark Matter | ``swift_type_dark_matter`` |
-+---------------------+------------------------+----------------------------+
-| ``/PartType2/`` | Ignored | |
-+---------------------+------------------------+----------------------------+
-| ``/PartType3/`` | Ignored | |
-+---------------------+------------------------+----------------------------+
-| ``/PartType4/`` | Stars | ``swift_type_star`` |
-+---------------------+------------------------+----------------------------+
-| ``/PartType5/`` | Black Holes | ``swift_type_black_hole`` |
-+---------------------+------------------------+----------------------------+
++---------------------+------------------------+----------------------------------------+
+| HDF5 Group Name | Physical Particle Type | In code ``enum part_type`` |
++=====================+========================+========================================+
+| ``/PartType0/`` | Gas | ``swift_type_gas`` |
++---------------------+------------------------+----------------------------------------+
+| ``/PartType1/`` | Dark Matter | ``swift_type_dark_matter`` |
++---------------------+------------------------+----------------------------------------+
+| ``/PartType2/`` | Background Dark Matter | ``swift_type_dark_matter_background`` |
++---------------------+------------------------+----------------------------------------+
+| ``/PartType3/`` | Ignored | |
++---------------------+------------------------+----------------------------------------+
+| ``/PartType4/`` | Stars | ``swift_type_star`` |
++---------------------+------------------------+----------------------------------------+
+| ``/PartType5/`` | Black Holes | ``swift_type_black_hole`` |
++---------------------+------------------------+----------------------------------------+
The last column in the table gives the ``enum`` value from ``part_type.h``
corresponding to a given entry in the files.
-Note that the only particles that have hydrodynamical forces calculated between
-them are those in ``PartType0``.
+Note that the only particles that have hydrodynamical forces calculated
+between them are those in ``PartType0``. The background dark matter
+particles are used for zoom-in simulations and can have different masses
+(and as a consequence softening length) within the ``/PartType2`` arrays.
Necessary Components
@@ -137,8 +139,8 @@ individual particle type (e.g. ``/PartType0/``) that have the following *dataset
within [0, L)^3 where L is the side-length of the simulation volume. In the
case of cosmological simulations, these are the co-moving positions.
+ ``Velocities``, an array of shape (N, 3) that is the cartesian velocities of
- the particles. When running cosmological simulations, these are the peculiar
- velocities. Note that this is different from GADGET which uses peculiar
+ the particles. When running cosmological simulations, these are the **peculiar
+ velocities**. Note that this is different from GADGET which uses peculiar
velocities divided by ``sqrt(a)`` (see below for a fix).
+ ``ParticleIDs``, an array of length N that are unique identifying numbers for
each particle. Note that these have to be unique to a particle, and cannot be
diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst
index b90c2c26b8fe01ad306106fd3fd2f65612002d8f..16938a0fa14d5cfa9b8bb14fa9d784eb0ed08c82 100644
--- a/doc/RTD/source/ParameterFiles/parameter_description.rst
+++ b/doc/RTD/source/ParameterFiles/parameter_description.rst
@@ -174,7 +174,7 @@ use the following parameters:
h: 0.6777
Omega_m: 0.307
Omega_lambda: 0.693
- Omega_b: 0.0455
+ Omega_b: 0.0482519
Omega_r: 0. # (Optional)
w_0: -1.0 # (Optional)
w_a: 0. # (Optional)
@@ -191,12 +191,32 @@ The behaviour of the self-gravity solver can be modified by the parameters
provided in the ``Gravity`` section. The theory document puts these parameters into the
context of the equations being solved. We give a brief overview here.
-* The Plummer-equivalent co-moving softening length used for all particles :math:`\epsilon_{com}`: ``comoving_softening``,
-* The Plummer-equivalent maximal physical softening length used for all particles :math:`\epsilon_{max}`: ``comoving_softening``,
-
-At any redshift :math:`z`, the Plummer-equivalent softening length used by the
-code will be :math:`\epsilon=\min(\epsilon_{max},
-\frac{\epsilon_{com}}{z+1})`. This is expressed in internal units.
+* The Plummer-equivalent co-moving softening length used for all dark matter particles :math:`\epsilon_{\rm com,DM}`: ``comoving_DM_softening``,
+* The Plummer-equivalent co-moving softening length used for all baryon particles (gas, stars, BHs) :math:`\epsilon_{\rm com,bar}`: ``comoving_baryon_softening``,
+* The Plummer-equivalent maximal physical softening length used for all dark matter particles :math:`\epsilon_{\rm max,DM}`: ``max_physical_DM_softening``,
+* The Plummer-equivalent maximal physical softening length used for all baryon particles (gas, stars, BHs) :math:`\epsilon_{\rm max,bar}`: ``max_physical_baryon_softening``,
+
+At any redshift :math:`z`, the Plummer-equivalent softening length used by
+the code will be :math:`\epsilon=\min(\epsilon_{max},
+\frac{\epsilon_{com}}{z+1})`. The same calculation is performed
+independently for the dark matter and baryon particles. All the softening
+quantities are expressed in internal units. Calculations that only involve
+DM or baryons can leave the unused quantities out of the parameter
+file. For non-cosmological runs, only the physical softening lengths need
+to be supplied.
+
+In case of zoom simulations, the softening of the additional, more massive, background
+particles is specified via the parameter
+``softening_ratio_background``. Since these particles will typically have
+different masses to degrade the resolution away from the zoom region, the
+particles won't have a single softening value. Instead, we specify the
+fraction of the mean inter-particle separation to use. The code will then
+derive the softening length of each particle assuming the mean density of
+the Universe. That is :math:`\epsilon_{\rm background} =
+f\sqrt[3]{\frac{m}{\Omega_m\rho_{\rm crit}}}`, where :math:`f` is the
+user-defined value (typically of order 0.05).
+
+The accuracy of the gravity calculation is governed by the following two parameters:
* The opening angle (multipole acceptance criterion) used in the FMM :math:`\theta`: ``theta``,
* The time-step size pre-factor :math:`\eta`: ``eta``,
@@ -236,15 +256,17 @@ simulation:
# Parameters for the self-gravity scheme for the EAGLE-100 box
Gravity:
- eta: 0.025
- theta: 0.7
- comoving_softening: 0.0026994 # 0.7 proper kpc at z=2.8.
- max_physical_softening: 0.0007 # 0.7 proper kpc
- rebuild_frequency: 0.01 # Default optional value
+ eta: 0.025
+ theta: 0.7
mesh_side_length: 512
- a_smooth: 1.25 # Default optional value
- r_cut_max: 4.5 # Default optional value
- r_cut_min: 0.1 # Default optional value
+ comoving_DM_softening: 0.0026994 # 0.7 proper kpc at z=2.8.
+ max_physical_DM_softening: 0.0007 # 0.7 proper kpc
+ comoving_baryon_softening: 0.0026994 # 0.7 proper kpc at z=2.8.
+ max_physical_baryon_softening: 0.0007 # 0.7 proper kpc
+ rebuild_frequency: 0.01 # Default optional value
+ a_smooth: 1.25 # Default optional value
+ r_cut_max: 4.5 # Default optional value
+ r_cut_min: 0.1 # Default optional value
.. _Parameters_SPH:
diff --git a/doc/RTD/source/Snapshots/index.rst b/doc/RTD/source/Snapshots/index.rst
index 32317aa19b8a216659317d4ae457a2dbfd040571..fa11058edfa202b548936dac3ca55ada86eb9e84 100644
--- a/doc/RTD/source/Snapshots/index.rst
+++ b/doc/RTD/source/Snapshots/index.rst
@@ -107,17 +107,19 @@ There are several groups that contain 'auxiliary' information, such as
the particles. The type use the naming convention of Gadget-2 (with
the OWLS and EAGLE extensions).
-+---------------------+------------------------+----------------------------+
-| HDF5 Group Name | Physical Particle Type | In code ``enum part_type`` |
-+=====================+========================+============================+
-| ``/PartType0/`` | Gas | ``swift_type_gas`` |
-+---------------------+------------------------+----------------------------+
-| ``/PartType1/`` | Dark Matter | ``swift_type_dark_matter`` |
-+---------------------+------------------------+----------------------------+
-| ``/PartType4/`` | Stars | ``swift_type_star`` |
-+---------------------+------------------------+----------------------------+
-| ``/PartType5/`` | Black Holes | ``swift_type_black_hole`` |
-+---------------------+------------------------+----------------------------+
++---------------------+------------------------+----------------------------------------+
+| HDF5 Group Name | Physical Particle Type | In code ``enum part_type`` |
++=====================+========================+========================================+
+| ``/PartType0/`` | Gas | ``swift_type_gas`` |
++---------------------+------------------------+----------------------------------------+
+| ``/PartType1/`` | Dark Matter | ``swift_type_dark_matter`` |
++---------------------+------------------------+----------------------------------------+
+| ``/PartType2/`` | Background Dark Matter | ``swift_type_dark_matter_background`` |
++---------------------+------------------------+----------------------------------------+
+| ``/PartType4/`` | Stars | ``swift_type_star`` |
++---------------------+------------------------+----------------------------------------+
+| ``/PartType5/`` | Black Holes | ``swift_type_black_hole`` |
++---------------------+------------------------+----------------------------------------+
The last column in the table gives the ``enum`` value from ``part_type.h``
corresponding to a given entry in the files.
diff --git a/examples/Cosmology/ConstantCosmoVolume/constant_volume.yml b/examples/Cosmology/ConstantCosmoVolume/constant_volume.yml
index ebfcc4ffd72121571fa1a69f900985917b440c65..a6ff72555ef68964508493856127d4cc739b7722 100644
--- a/examples/Cosmology/ConstantCosmoVolume/constant_volume.yml
+++ b/examples/Cosmology/ConstantCosmoVolume/constant_volume.yml
@@ -49,6 +49,8 @@ Gravity:
mesh_side_length: 32
eta: 0.025
theta: 0.3
- comoving_softening: 0.08 # 80 kpc = 1/25 of mean inter-particle separation
- max_physical_softening: 0.08 # 80 kpc = 1/25 of mean inter-particle separation
+ comoving_DM_softening: 0.08 # 80 kpc = 1/25 of mean inter-particle separation
+ max_physical_DM_softening: 0.08 # 80 kpc = 1/25 of mean inter-particle separation
+ comoving_baryon_softening: 0.08 # 80 kpc = 1/25 of mean inter-particle separation
+ max_physical_baryon_softening: 0.08 # 80 kpc = 1/25 of mean inter-particle separation
diff --git a/examples/Cosmology/ZeldovichPancake_3D/zeldovichPancake.yml b/examples/Cosmology/ZeldovichPancake_3D/zeldovichPancake.yml
index 6a7c5166635b7fa0ed5f69c41461d867c3b254ad..d43c78972b0bc8d1f250b95190dafef305abca3f 100644
--- a/examples/Cosmology/ZeldovichPancake_3D/zeldovichPancake.yml
+++ b/examples/Cosmology/ZeldovichPancake_3D/zeldovichPancake.yml
@@ -51,5 +51,7 @@ Gravity:
eta: 0.025
theta: 0.3
r_cut_max: 5.
- comoving_softening: 0.001
- max_physical_softening: 0.001
+ comoving_DM_softening: 0.001
+ max_physical_DM_softening: 0.001
+ comoving_baryon_softening: 0.001
+ max_physical_baryon_softening: 0.001
diff --git a/examples/EAGLE_DMO_low_z/EAGLE_DMO_100/eagle_100.yml b/examples/EAGLE_DMO_low_z/EAGLE_DMO_100/eagle_100.yml
index 5a3066195647b79eeb6a6d67d037d15ce8370c39..5d449cfceb589e7580ba5f52375d115678f51321 100644
--- a/examples/EAGLE_DMO_low_z/EAGLE_DMO_100/eagle_100.yml
+++ b/examples/EAGLE_DMO_low_z/EAGLE_DMO_100/eagle_100.yml
@@ -23,7 +23,7 @@ TimeIntegration:
dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
Scheduler:
- max_top_level_cells: 80
+ max_top_level_cells: 64
# Parameters governing the snapshots
Snapshots:
@@ -43,8 +43,9 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.85 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
+ mesh_side_length: 512
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
# Parameters related to the initial conditions
InitialConditions:
diff --git a/examples/EAGLE_DMO_low_z/EAGLE_DMO_12/eagle_12.yml b/examples/EAGLE_DMO_low_z/EAGLE_DMO_12/eagle_12.yml
index 0660d98e87adfae62a2d795efec7ad6509cc1354..f232bb7e930a4977b65686424d04a7eed254761d 100644
--- a/examples/EAGLE_DMO_low_z/EAGLE_DMO_12/eagle_12.yml
+++ b/examples/EAGLE_DMO_low_z/EAGLE_DMO_12/eagle_12.yml
@@ -41,12 +41,11 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
- mesh_side_length: 32
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
mesh_side_length: 32
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
# Parameters related to the initial conditions
InitialConditions:
diff --git a/examples/EAGLE_DMO_low_z/EAGLE_DMO_25/eagle_25.yml b/examples/EAGLE_DMO_low_z/EAGLE_DMO_25/eagle_25.yml
index 558c68ffaad204ebbe1d5781f945f0d95108d227..20f4c1eef17fafd8a87d1d251147438fc201a355 100644
--- a/examples/EAGLE_DMO_low_z/EAGLE_DMO_25/eagle_25.yml
+++ b/examples/EAGLE_DMO_low_z/EAGLE_DMO_25/eagle_25.yml
@@ -43,9 +43,9 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
- mesh_side_length: 32
+ mesh_side_length: 64
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
# Parameters related to the initial conditions
InitialConditions:
diff --git a/examples/EAGLE_DMO_low_z/EAGLE_DMO_50/eagle_50.yml b/examples/EAGLE_DMO_low_z/EAGLE_DMO_50/eagle_50.yml
index 3cab2b1dc869b5187cf647caa7893281b783591a..0ccf061c2e942126bbc0c2c69998bbca6a9b824f 100644
--- a/examples/EAGLE_DMO_low_z/EAGLE_DMO_50/eagle_50.yml
+++ b/examples/EAGLE_DMO_low_z/EAGLE_DMO_50/eagle_50.yml
@@ -23,7 +23,7 @@ TimeIntegration:
dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
Scheduler:
- max_top_level_cells: 20
+ max_top_level_cells: 32
# Parameters governing the snapshots
Snapshots:
@@ -42,9 +42,9 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
- mesh_side_length: 64
+ mesh_side_length: 128
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
# Parameters related to the initial conditions
InitialConditions:
diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
index e94dce66f80b9b7c669e9095b7ce6f76aaf84551..93c8c740f58efb23a017a7d229f81a685e837b1a 100644
--- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
@@ -39,9 +39,11 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
mesh_side_length: 64
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
+ comoving_baryon_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_baryon_softening: 0.0007 # Max physical DM softening length (in internal units).
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
index 4e6f7f6c49040fb2f39f443c7c1c7235570566f3..d5307533c0ffccae7644f06dbe33c27bf46f4114 100644
--- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
@@ -39,9 +39,11 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
mesh_side_length: 128
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
+ comoving_baryon_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_baryon_softening: 0.0007 # Max physical DM softening length (in internal units).
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
index 9317e9603b106141e9c7a82b9cd2fc36f63698a8..f757dc9dcf104237c6ecc5e472d29f79375a1d53 100644
--- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
@@ -39,9 +39,11 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
mesh_side_length: 256
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
+ comoving_baryon_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_baryon_softening: 0.0007 # Max physical DM softening length (in internal units).
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_low_z/EAGLE_100/eagle_100.yml b/examples/EAGLE_low_z/EAGLE_100/eagle_100.yml
index 746362eb7d573eed75018423b8e33f6835916ccd..f77036b5d55f33b4fd3f42c7bea0ccc124003a40 100644
--- a/examples/EAGLE_low_z/EAGLE_100/eagle_100.yml
+++ b/examples/EAGLE_low_z/EAGLE_100/eagle_100.yml
@@ -1,3 +1,6 @@
+MetaData:
+ run_name: EAGLE-L0100N1504-Ref
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
@@ -42,9 +45,11 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.85 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
mesh_side_length: 256
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
+ comoving_baryon_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_baryon_softening: 0.0007 # Max physical DM softening length (in internal units).
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index da6dfc0d094e45146e8bed3e536ca7a0874093fb..f470c691a5a76207998f6d854f6e8d44f0a1aebb 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -1,3 +1,6 @@
+MetaData:
+ run_name: EAGLE-L0012N0188-Ref
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
@@ -43,9 +46,11 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
mesh_side_length: 32
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
+ comoving_baryon_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_baryon_softening: 0.0007 # Max physical DM softening length (in internal units).
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index 1c2f62ef75209d8e6b3ed787b97d458e6b567dc6..07685bf783b34b2872df4a32610fa791db01cded 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -1,3 +1,6 @@
+MetaData:
+ run_name: EAGLE-L0025N0376-Ref
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
@@ -49,10 +52,13 @@ 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: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
mesh_side_length: 64
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
+ comoving_baryon_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_baryon_softening: 0.0007 # Max physical DM softening length (in internal units).
+
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index 14844f49a5fe56607e0bd95dd57d44c42cb43f9d..943c64c7a29fd87b9b5f78a4edded6b14e0f3c57 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -1,3 +1,6 @@
+MetaData:
+ run_name: EAGLE-L0050N0752-Ref
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
@@ -42,9 +45,11 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
mesh_side_length: 128
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
+ comoving_baryon_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_baryon_softening: 0.0007 # Max physical DM softening length (in internal units).
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index 9aef75d741efa148f2397ba172d1ca6c86dd20d9..48a825750fd2a927ba08dfc5a8a4607a490fe0d8 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -1,3 +1,6 @@
+MetaData:
+ run_name: EAGLE-L0006N0094-Ref
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
@@ -52,9 +55,11 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
mesh_side_length: 16
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
+ comoving_baryon_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_baryon_softening: 0.0007 # Max physical DM softening length (in internal units).
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/GiantImpacts/GiantImpact/uranus_1e6.yml b/examples/GiantImpacts/GiantImpact/uranus_1e6.yml
index 355748d847097623f171078c2ca8372e06a06efa..b98467616f659babb488352c87697a43fedcc0db 100644
--- a/examples/GiantImpacts/GiantImpact/uranus_1e6.yml
+++ b/examples/GiantImpacts/GiantImpact/uranus_1e6.yml
@@ -48,8 +48,8 @@ SPH:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.003 # Comoving softening length (in internal units).
- max_physical_softening: 0.003 # Physical softening length (in internal units).
+ comoving_baryon_softening: 0.003 # Comoving softening length (in internal units).
+ max_physical_baryon_softening: 0.003 # Physical softening length (in internal units).
# Parameters for the task scheduling
Scheduler:
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_dmparticles/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_dmparticles/isolated_galaxy.yml
index dccfb28a3f1c888d2a83b5e28b759a30a6928754..27ab01d984319fea68d4d1ae8fb435a6adf895ce 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_dmparticles/isolated_galaxy.yml
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_dmparticles/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.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
+ max_physical_DM_softening: 0.7 # 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:
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
index 4d3485d838f56067cbb8493c79d06dff8a3dceba..79fe5682692347127081f021ed6930df57bbfa02 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
+++ b/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:
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_potential/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_potential/isolated_galaxy.yml
index deee132ee38ae5e04397839a21a677f4851e6bac..3bd743c2ec329057c7f63f987e4809744a07f7ba 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_potential/isolated_galaxy.yml
+++ b/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:
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml
index d347e1ad5bd8ddf6c132fade07e5d50688a451b0..d917f926724c022cd15524058ddde2a7466acaab 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml
+++ b/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:
diff --git a/examples/PMillennium/PMillennium-1536/p-mill-1536.yml b/examples/PMillennium/PMillennium-1536/p-mill-1536.yml
index 770cc4a03e8912eed51d1851b81578496a42070c..f343650452a24f620edef88e42666e09213dbd64 100644
--- a/examples/PMillennium/PMillennium-1536/p-mill-1536.yml
+++ b/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:
diff --git a/examples/PMillennium/PMillennium-384/p-mill-384.yml b/examples/PMillennium/PMillennium-384/p-mill-384.yml
index 428f893300bf3431f20518a6a1259651b3c67226..0e68969d0b590cec8058805e4644b2763c543be1 100644
--- a/examples/PMillennium/PMillennium-384/p-mill-384.yml
+++ b/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:
diff --git a/examples/PMillennium/PMillennium-768/p-mill-768.yml b/examples/PMillennium/PMillennium-768/p-mill-768.yml
index 339fad2c7059d8c8bcab2878fc8958a3bd1977d7..1cd9e63b1f03ac65baf72701de276b8ff43c9575 100644
--- a/examples/PMillennium/PMillennium-768/p-mill-768.yml
+++ b/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:
diff --git a/examples/SantaBarbara/SantaBarbara-128/santa_barbara.yml b/examples/SantaBarbara/SantaBarbara-128/santa_barbara.yml
index 1221cefb92e51e54bfe73a0414b95a89fcda592a..e57617e149cb9f302d8fe08bd93d660059caa0ef 100644
--- a/examples/SantaBarbara/SantaBarbara-128/santa_barbara.yml
+++ b/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:
diff --git a/examples/SantaBarbara/SantaBarbara-256/santa_barbara.yml b/examples/SantaBarbara/SantaBarbara-256/santa_barbara.yml
index 04bb18191b750aa996e43b58b24451135bceedca..c95914d58ef01ac48f5db83ea4b3f3a1f93404f9 100644
--- a/examples/SantaBarbara/SantaBarbara-256/santa_barbara.yml
+++ b/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:
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_DM/small_cosmo_volume_dm.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_DM/small_cosmo_volume_dm.yml
index ebe3a78ee0d03eb53752b1dfa8fa749931a754a9..975ceac9c3c9404dc2c194abc4de40bf7d529a11 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_DM/small_cosmo_volume_dm.yml
+++ b/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
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_VELOCIraptor/small_cosmo_volume.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_VELOCIraptor/small_cosmo_volume.yml
index 15007ca8d39a328166a208a2d4cbbc6aea580009..f81c38a85f04b074f4cf5598bc3d09716b49f790 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_VELOCIraptor/small_cosmo_volume.yml
+++ b/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
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
index 5050f9e05e278df25b050c10d9897e29d4e21804..aea198bd789b48d608c6fc3fdc493303d6d9afd2 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
+++ b/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
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/small_cosmo_volume.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/small_cosmo_volume.yml
index e0d633079e941ade161b7e2fde0fbc063cbac254..5904a412f7bf63336de976baa6da4506bae4e36c 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/small_cosmo_volume.yml
+++ b/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
diff --git a/examples/main.c b/examples/main.c
index 491b1ac1e72895ba7f9028d7c0b1bf9181f70360..62e5e85ba13527f1119b73234e2ceb1122f9f397 100644
--- a/examples/main.c
+++ b/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)...",
diff --git a/examples/main_fof.c b/examples/main_fof.c
index 1672cdacd5322d8c18922c1b9784f779e776c0ac..c96cb1267e8c3ab6890de121a9d7fb3fb16443df 100644
--- a/examples/main_fof.c
+++ b/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)...",
diff --git a/examples/nIFTyCluster/Baryonic/clean.sh b/examples/nIFTyCluster/Baryonic/clean.sh
old mode 100644
new mode 100755
diff --git a/examples/nIFTyCluster/Baryonic/create_ics.sh b/examples/nIFTyCluster/Baryonic/create_ics.sh
old mode 100644
new mode 100755
diff --git a/examples/nIFTyCluster/Baryonic/getIC.sh b/examples/nIFTyCluster/Baryonic/getIC.sh
old mode 100644
new mode 100755
diff --git a/examples/nIFTyCluster/Baryonic/nifty.yml b/examples/nIFTyCluster/Baryonic/nifty.yml
index d79d6853a59722a19b740b125d21883c2c06aa2b..d520d5bfba6ca466c30a0bd234a5457d7edd8fdf 100644
--- a/examples/nIFTyCluster/Baryonic/nifty.yml
+++ b/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:
diff --git a/examples/nIFTyCluster/Baryonic/run.sh b/examples/nIFTyCluster/Baryonic/run.sh
old mode 100644
new mode 100755
diff --git a/examples/nIFTyCluster/Baryonic/setup.sh b/examples/nIFTyCluster/Baryonic/setup.sh
old mode 100644
new mode 100755
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index aabf9ac1078ad41288f1011d7bfe94626deefce3..7d8a5886ebcde64a39defa24b9f1e9fbd2232642 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -21,7 +21,7 @@ Cosmology:
a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
Omega_r: 0. # (Optional) Radiation density parameter
w_0: -1.0 # (Optional) Dark-energy equation-of-state parameter at z=0.
w_a: 0. # (Optional) Dark-energy equation-of-state time evolution parameter.
@@ -58,15 +58,18 @@ Stars:
# 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.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
- rebuild_frequency: 0.01 # (Optional) Frequency of the gravity-tree rebuild in units of the number of g-particles (this is the default value).
- a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
- r_cut_max: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
- r_cut_min: 0.1 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
+ 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_DM_softening: 0.0026994 # Comoving Plummer-equivalent softening length for DM particles (in internal units).
+ max_physical_DM_softening: 0.0007 # Maximal Plummer-equivalent softening length in physical coordinates for DM particles (in internal units).
+ comoving_baryon_softening: 0.0026994 # Comoving Plummer-equivalent softening length for baryon particles (in internal units).
+ max_physical_baryon_softening: 0.0007 # Maximal Plummer-equivalent softening length in physical coordinates for baryon particles (in internal units).
+ softening_ratio_background: 0.04 # Fraction of the mean inter-particle separation to use as Plummer-equivalent softening for the background DM particles.
+ rebuild_frequency: 0.01 # (Optional) Frequency of the gravity-tree rebuild in units of the number of g-particles (this is the default value).
+ a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
+ r_cut_max: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
+ r_cut_min: 0.1 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
# Parameters for the Friends-Of-Friends algorithm
FOF:
diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h
index 7e9852dee582ab49643368c0a7b0d79ecd6b5047..3e2238dc2a05f73dc9020b66f754e26a5976572a 100644
--- a/src/black_holes/EAGLE/black_holes_iact.h
+++ b/src/black_holes/EAGLE/black_holes_iact.h
@@ -145,7 +145,7 @@ runner_iact_nonsym_bh_gas_swallow(
const float max_dist_repos2 =
kernel_gravity_softening_plummer_equivalent_inv *
kernel_gravity_softening_plummer_equivalent_inv * 9.f *
- grav_props->epsilon_cur2;
+ grav_props->epsilon_baryon_cur * grav_props->epsilon_baryon_cur;
/* This gas neighbour is close enough that we can consider it's potential
for repositioning */
@@ -250,7 +250,7 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
const float max_dist_repos2 =
kernel_gravity_softening_plummer_equivalent_inv *
kernel_gravity_softening_plummer_equivalent_inv * 9.f *
- grav_props->epsilon_cur2;
+ grav_props->epsilon_baryon_cur * grav_props->epsilon_baryon_cur;
/* This gas neighbour is close enough that we can consider it's potential
for repositioning */
@@ -285,7 +285,7 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
const float max_dist_merge2 =
kernel_gravity_softening_plummer_equivalent_inv *
kernel_gravity_softening_plummer_equivalent_inv * 9.f *
- grav_props->epsilon_cur2;
+ grav_props->epsilon_baryon_cur * grav_props->epsilon_baryon_cur;
const float G_Newton = grav_props->G_Newton;
diff --git a/src/black_holes/EAGLE/black_holes_part.h b/src/black_holes/EAGLE/black_holes_part.h
index cbfa0fcaf560d1c771d8595ec84ad08e9ba3608e..4d8a9c8e177cc64d7e3191c046ff95c374b2d4a4 100644
--- a/src/black_holes/EAGLE/black_holes_part.h
+++ b/src/black_holes/EAGLE/black_holes_part.h
@@ -19,6 +19,7 @@
#ifndef SWIFT_EAGLE_BLACK_HOLE_PART_H
#define SWIFT_EAGLE_BLACK_HOLE_PART_H
+#include "black_holes_struct.h"
#include "chemistry_struct.h"
#include "timeline.h"
diff --git a/src/black_holes/EAGLE/black_holes_struct.h b/src/black_holes/EAGLE/black_holes_struct.h
index 98f5380e9f32ced36bcf2e4e3f02eaeb89db7218..177a54b921bec24f357756e4f777b0ecf1781c30 100644
--- a/src/black_holes/EAGLE/black_holes_struct.h
+++ b/src/black_holes/EAGLE/black_holes_struct.h
@@ -19,6 +19,8 @@
#ifndef SWIFT_BLACK_HOLES_STRUCT_EAGLE_H
#define SWIFT_BLACK_HOLES_STRUCT_EAGLE_H
+#include "inline.h"
+
/**
* @brief Black holes-related fields carried by each *gas* particle.
*/
diff --git a/src/cell.c b/src/cell.c
index f2b38ee1a7d084a17a8e296f18136309136564a8..4b9746e92e31f2a5ae6c4b5d01ade8d83939e412 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -2178,8 +2178,11 @@ void cell_reset_task_counters(struct cell *c) {
*
* @param c The #cell.
* @param ti_current The current integer time.
+ * @param grav_props The properties of the gravity scheme.
*/
-void cell_make_multipoles(struct cell *c, integertime_t ti_current) {
+void cell_make_multipoles(struct cell *c, integertime_t ti_current,
+ const struct gravity_props *const grav_props) {
+
/* Reset everything */
gravity_reset(c->grav.multipole);
@@ -2187,7 +2190,7 @@ void cell_make_multipoles(struct cell *c, integertime_t ti_current) {
/* Start by recursing */
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL)
- cell_make_multipoles(c->progeny[k], ti_current);
+ cell_make_multipoles(c->progeny[k], ti_current, grav_props);
}
/* Compute CoM of all progenies */
@@ -2248,7 +2251,7 @@ void cell_make_multipoles(struct cell *c, integertime_t ti_current) {
} else {
if (c->grav.count > 0) {
- gravity_P2M(c->grav.multipole, c->grav.parts, c->grav.count);
+ gravity_P2M(c->grav.multipole, c->grav.parts, c->grav.count, grav_props);
const double dx =
c->grav.multipole->CoM[0] > c->loc[0] + c->width[0] * 0.5
? c->grav.multipole->CoM[0] - c->loc[0]
@@ -2324,8 +2327,11 @@ void cell_check_foreign_multipole(const struct cell *c) {
* recursively computed one.
*
* @param c Cell to act upon
+ * @param grav_props The properties of the gravity scheme.
*/
-void cell_check_multipole(struct cell *c) {
+void cell_check_multipole(struct cell *c,
+ const struct gravity_props *const grav_props) {
+
#ifdef SWIFT_DEBUG_CHECKS
struct gravity_tensors ma;
const double tolerance = 1e-3; /* Relative */
@@ -2333,11 +2339,12 @@ void cell_check_multipole(struct cell *c) {
/* First recurse */
if (c->split)
for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL) cell_check_multipole(c->progeny[k]);
+ if (c->progeny[k] != NULL)
+ cell_check_multipole(c->progeny[k], grav_props);
if (c->grav.count > 0) {
/* Brute-force calculation */
- gravity_P2M(&ma, c->grav.parts, c->grav.count);
+ gravity_P2M(&ma, c->grav.parts, c->grav.count, grav_props);
/* Now compare the multipole expansion */
if (!gravity_multipole_equal(&ma, c->grav.multipole, tolerance)) {
@@ -5450,6 +5457,9 @@ void cell_remove_gpart(const struct engine *e, struct cell *c,
if (c->nodeID != e->nodeID)
error("Can't remove a particle in a foreign cell.");
+ if (gp->type == swift_type_dark_matter_background)
+ error("Can't remove a DM background particle!");
+
/* Mark the particle as inhibited */
gp->time_bin = time_bin_inhibited;
@@ -5900,7 +5910,11 @@ int cell_can_use_pair_mm(const struct cell *ci, const struct cell *cj,
}
const double r2 = dx * dx + dy * dy + dz * dz;
- return gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2);
+ const double epsilon_i = multi_i->m_pole.max_softening;
+ const double epsilon_j = multi_j->m_pole.max_softening;
+
+ return gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2,
+ epsilon_i, epsilon_j);
}
/**
@@ -5962,6 +5976,9 @@ int cell_can_use_pair_mm_rebuild(const struct cell *ci, const struct cell *cj,
}
const double r2 = dx * dx + dy * dy + dz * dz;
+ const double epsilon_i = multi_i->m_pole.max_softening;
+ const double epsilon_j = multi_j->m_pole.max_softening;
+
return gravity_M2L_accept(multi_i->r_max_rebuild, multi_j->r_max_rebuild,
- theta_crit2, r2);
+ theta_crit2, r2, epsilon_i, epsilon_j);
}
diff --git a/src/cell.h b/src/cell.h
index 830b281b5a3766f4842ea8cd8276f41a8ef00b66..8067a3189818ab8738de848ea698fbf25c78ebba 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -854,8 +854,10 @@ int cell_link_foreign_gparts(struct cell *c, struct gpart *gparts);
int cell_count_parts_for_tasks(const struct cell *c);
int cell_count_gparts_for_tasks(const struct cell *c);
void cell_clean_links(struct cell *c, void *data);
-void cell_make_multipoles(struct cell *c, integertime_t ti_current);
-void cell_check_multipole(struct cell *c);
+void cell_make_multipoles(struct cell *c, integertime_t ti_current,
+ const struct gravity_props *const grav_props);
+void cell_check_multipole(struct cell *c,
+ const struct gravity_props *const grav_props);
void cell_check_foreign_multipole(const struct cell *c);
void cell_clean(struct cell *c);
void cell_check_part_drift_point(struct cell *c, void *data);
diff --git a/src/common_io.c b/src/common_io.c
index 945badd14cb814e5c06c5e0edf2b983669259850..6f090ee01d1c870408223f90dba0770b5f76f567 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -24,7 +24,7 @@
/* This object's header. */
#include "common_io.h"
-/* First common header */
+/* Pre-inclusion as needed in other headers */
#include "engine.h"
/* Local includes. */
@@ -422,6 +422,21 @@ static long long cell_count_non_inhibited_dark_matter(const struct cell* c) {
return count;
}
+static long long cell_count_non_inhibited_background_dark_matter(
+ const struct cell* c) {
+ const int total_count = c->grav.count;
+ struct gpart* gparts = c->grav.parts;
+ long long count = 0;
+ for (int i = 0; i < total_count; ++i) {
+ if ((gparts[i].time_bin != time_bin_inhibited) &&
+ (gparts[i].time_bin != time_bin_not_created) &&
+ (gparts[i].type == swift_type_dark_matter_background)) {
+ ++count;
+ }
+ }
+ return count;
+}
+
static long long cell_count_non_inhibited_stars(const struct cell* c) {
const int total_count = c->stars.count;
struct spart* sparts = c->stars.parts;
@@ -463,30 +478,36 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
centres = (double*)malloc(3 * nr_cells * sizeof(double));
/* Count of particles in each cell */
- long long *count_part = NULL, *count_gpart = NULL, *count_spart = NULL,
+ long long *count_part = NULL, *count_gpart = NULL,
+ *count_background_gpart = NULL, *count_spart = NULL,
*count_bpart = NULL;
count_part = (long long*)malloc(nr_cells * sizeof(long long));
count_gpart = (long long*)malloc(nr_cells * sizeof(long long));
+ count_background_gpart = (long long*)malloc(nr_cells * sizeof(long long));
count_spart = (long long*)malloc(nr_cells * sizeof(long long));
count_bpart = (long long*)malloc(nr_cells * sizeof(long long));
/* Global offsets of particles in each cell */
- long long *offset_part = NULL, *offset_gpart = NULL, *offset_spart = NULL,
+ long long *offset_part = NULL, *offset_gpart = NULL,
+ *offset_background_gpart = NULL, *offset_spart = NULL,
*offset_bpart = NULL;
offset_part = (long long*)malloc(nr_cells * sizeof(long long));
offset_gpart = (long long*)malloc(nr_cells * sizeof(long long));
+ offset_background_gpart = (long long*)malloc(nr_cells * sizeof(long long));
offset_spart = (long long*)malloc(nr_cells * sizeof(long long));
offset_bpart = (long long*)malloc(nr_cells * sizeof(long long));
/* Offsets of the 0^th element */
offset_part[0] = 0;
offset_gpart[0] = 0;
+ offset_background_gpart[0] = 0;
offset_spart[0] = 0;
offset_bpart[0] = 0;
/* Collect the cell information of *local* cells */
long long local_offset_part = 0;
long long local_offset_gpart = 0;
+ long long local_offset_background_gpart = 0;
long long local_offset_spart = 0;
long long local_offset_bpart = 0;
for (int i = 0; i < nr_cells; ++i) {
@@ -501,6 +522,8 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
/* Count real particles that will be written */
count_part[i] = cell_count_non_inhibited_gas(&cells_top[i]);
count_gpart[i] = cell_count_non_inhibited_dark_matter(&cells_top[i]);
+ count_background_gpart[i] =
+ cell_count_non_inhibited_background_dark_matter(&cells_top[i]);
count_spart[i] = cell_count_non_inhibited_stars(&cells_top[i]);
count_bpart[i] = cell_count_non_inhibited_black_holes(&cells_top[i]);
@@ -509,12 +532,16 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
offset_part[i] = local_offset_part + global_offsets[swift_type_gas];
offset_gpart[i] =
local_offset_gpart + global_offsets[swift_type_dark_matter];
+ offset_background_gpart[i] =
+ local_offset_background_gpart +
+ global_offsets[swift_type_dark_matter_background];
offset_spart[i] = local_offset_spart + global_offsets[swift_type_stars];
offset_bpart[i] =
local_offset_bpart + global_offsets[swift_type_black_hole];
local_offset_part += count_part[i];
local_offset_gpart += count_gpart[i];
+ local_offset_background_gpart += count_background_gpart[i];
local_offset_spart += count_spart[i];
local_offset_bpart += count_bpart[i];
@@ -528,11 +555,13 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
count_part[i] = 0;
count_gpart[i] = 0;
+ count_background_gpart[i] = 0;
count_spart[i] = 0;
count_bpart[i] = 0;
offset_part[i] = 0;
offset_gpart[i] = 0;
+ offset_background_gpart[i] = 0;
offset_spart[i] = 0;
offset_bpart[i] = 0;
}
@@ -548,7 +577,6 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
MPI_Reduce(count_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
MPI_COMM_WORLD);
}
-
if (nodeID == 0) {
MPI_Reduce(MPI_IN_PLACE, count_gpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
0, MPI_COMM_WORLD);
@@ -556,6 +584,13 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
MPI_Reduce(count_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
MPI_COMM_WORLD);
}
+ if (nodeID == 0) {
+ MPI_Reduce(MPI_IN_PLACE, count_background_gpart, nr_cells,
+ MPI_LONG_LONG_INT, MPI_BOR, 0, MPI_COMM_WORLD);
+ } else {
+ MPI_Reduce(count_background_gpart, NULL, nr_cells, MPI_LONG_LONG_INT,
+ MPI_BOR, 0, MPI_COMM_WORLD);
+ }
if (nodeID == 0) {
MPI_Reduce(MPI_IN_PLACE, count_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
0, MPI_COMM_WORLD);
@@ -585,6 +620,13 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
MPI_Reduce(offset_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
MPI_COMM_WORLD);
}
+ if (nodeID == 0) {
+ MPI_Reduce(MPI_IN_PLACE, offset_background_gpart, nr_cells,
+ MPI_LONG_LONG_INT, MPI_BOR, 0, MPI_COMM_WORLD);
+ } else {
+ MPI_Reduce(offset_background_gpart, NULL, nr_cells, MPI_LONG_LONG_INT,
+ MPI_BOR, 0, MPI_COMM_WORLD);
+ }
if (nodeID == 0) {
MPI_Reduce(MPI_IN_PLACE, offset_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
0, MPI_COMM_WORLD);
@@ -700,6 +742,28 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
H5Sclose(h_space);
}
+ if (global_counts[swift_type_dark_matter_background] > 0) {
+
+ shape[0] = nr_cells;
+ shape[1] = 1;
+ h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0)
+ error("Error while creating data space for background DM offsets");
+ h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
+ if (h_err < 0)
+ error(
+ "Error while changing shape of background DM offsets data space.");
+ h_data = H5Dcreate(h_subgrp, "PartType2", io_hdf5_type(LONGLONG), h_space,
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_data < 0)
+ error("Error while creating dataspace for background DM offsets.");
+ h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
+ H5P_DEFAULT, offset_background_gpart);
+ if (h_err < 0) error("Error while writing background DM offsets.");
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+ }
+
if (global_counts[swift_type_stars] > 0) {
shape[0] = nr_cells;
@@ -786,6 +850,27 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
H5Sclose(h_space);
}
+ if (global_counts[swift_type_dark_matter_background] > 0) {
+
+ shape[0] = nr_cells;
+ shape[1] = 1;
+ h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0)
+ error("Error while creating data space for background DM counts");
+ h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
+ if (h_err < 0)
+ error("Error while changing shape of background DM counts data space.");
+ h_data = H5Dcreate(h_subgrp, "PartType2", io_hdf5_type(LONGLONG), h_space,
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_data < 0)
+ error("Error while creating dataspace for background DM counts.");
+ h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
+ H5P_DEFAULT, count_background_gpart);
+ if (h_err < 0) error("Error while writing background DM counts.");
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+ }
+
if (global_counts[swift_type_stars] > 0) {
shape[0] = nr_cells;
@@ -834,10 +919,12 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
free(centres);
free(count_part);
free(count_gpart);
+ free(count_background_gpart);
free(count_spart);
free(count_bpart);
free(offset_part);
free(offset_gpart);
+ free(offset_background_gpart);
free(offset_spart);
free(offset_bpart);
}
@@ -1313,29 +1400,6 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
threadpool_map((struct threadpool*)&e->threadpool,
io_convert_gpart_f_mapper, temp_f, N, copySize, 0,
(void*)&props);
- } else if (props.convert_gpart_i != NULL) {
-
- /* Prepare some parameters */
- int* temp_i = (int*)temp;
- props.start_temp_i = (int*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_gpart_i_mapper, temp_i, N, copySize, 0,
- (void*)&props);
-
- } else if (props.convert_gpart_i != NULL) {
-
- /* Prepare some parameters */
- int* temp_i = (int*)temp;
- props.start_temp_i = (int*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_gpart_i_mapper, temp_i, N, copySize, 0,
- (void*)&props);
} else if (props.convert_gpart_i != NULL) {
@@ -1384,29 +1448,6 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
threadpool_map((struct threadpool*)&e->threadpool,
io_convert_spart_f_mapper, temp_f, N, copySize, 0,
(void*)&props);
- } else if (props.convert_spart_i != NULL) {
-
- /* Prepare some parameters */
- int* temp_i = (int*)temp;
- props.start_temp_i = (int*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_spart_i_mapper, temp_i, N, copySize, 0,
- (void*)&props);
-
- } else if (props.convert_spart_i != NULL) {
-
- /* Prepare some parameters */
- int* temp_i = (int*)temp;
- props.start_temp_i = (int*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_spart_i_mapper, temp_i, N, copySize, 0,
- (void*)&props);
} else if (props.convert_spart_i != NULL) {
@@ -1551,6 +1592,56 @@ void io_prepare_dm_gparts(struct threadpool* tp, struct gpart* const gparts,
sizeof(struct gpart), 0, NULL);
}
+void io_prepare_dm_background_gparts_mapper(void* restrict data, int Ndm,
+ void* dummy) {
+
+ struct gpart* restrict gparts = (struct gpart*)data;
+
+ /* Let's give all these gparts a negative id */
+ for (int i = 0; i < Ndm; ++i) {
+
+ /* Negative ids are not allowed */
+ if (gparts[i].id_or_neg_offset < 0)
+ error("Negative ID for DM particle %i: ID=%lld", i,
+ gparts[i].id_or_neg_offset);
+
+ /* Set gpart type */
+ gparts[i].type = swift_type_dark_matter_background;
+ }
+}
+
+/**
+ * @brief Prepare the DM backgorund particles (in gparts) read in
+ * for the addition of the other particle types
+ *
+ * This function assumes that the DM particles are all at the start of the
+ * gparts array and that the background particles directly follow them.
+ *
+ * @param tp The current #threadpool.
+ * @param gparts The array of #gpart freshly read in.
+ * @param Ndm The number of DM particles read in.
+ */
+void io_prepare_dm_background_gparts(struct threadpool* tp,
+ struct gpart* const gparts, size_t Ndm) {
+
+ threadpool_map(tp, io_prepare_dm_background_gparts_mapper, gparts, Ndm,
+ sizeof(struct gpart), 0, NULL);
+}
+
+size_t io_count_dm_background_gparts(const struct gpart* const gparts,
+ const size_t Ndm) {
+
+ swift_declare_aligned_ptr(const struct gpart, gparts_array, gparts,
+ gpart_align);
+
+ size_t count = 0;
+ for (size_t i = 0; i < Ndm; ++i) {
+ if (gparts[i].type == swift_type_dark_matter_background) ++count;
+ }
+
+ return count;
+}
+
struct duplication_data {
struct part* parts;
@@ -1846,7 +1937,8 @@ void io_collect_bparts_to_write(const struct bpart* restrict bparts,
}
/**
- * @brief Copy every non-inhibited DM #gpart into the gparts_written array.
+ * @brief Copy every non-inhibited regulat DM #gpart into the gparts_written
+ * array.
*
* @param gparts The array of #gpart containing all particles.
* @param vr_data The array of gpart-related VELOCIraptor output.
@@ -1888,6 +1980,50 @@ void io_collect_gparts_to_write(
count, Ngparts_written);
}
+/**
+ * @brief Copy every non-inhibited background DM #gpart into the gparts_written
+ * array.
+ *
+ * @param gparts The array of #gpart containing all particles.
+ * @param vr_data The array of gpart-related VELOCIraptor output.
+ * @param gparts_written The array of #gpart to fill with particles we want to
+ * write.
+ * @param vr_data_written The array of gpart-related VELOCIraptor with particles
+ * we want to write.
+ * @param Ngparts The total number of #part.
+ * @param Ngparts_written The total number of #part to write.
+ * @param with_stf Are we running with STF? i.e. do we want to collect vr data?
+ */
+void io_collect_gparts_background_to_write(
+ const struct gpart* restrict gparts,
+ const struct velociraptor_gpart_data* restrict vr_data,
+ struct gpart* restrict gparts_written,
+ struct velociraptor_gpart_data* restrict vr_data_written,
+ const size_t Ngparts, const size_t Ngparts_written, const int with_stf) {
+
+ size_t count = 0;
+
+ /* Loop over all parts */
+ for (size_t i = 0; i < Ngparts; ++i) {
+
+ /* And collect the ones that have not been removed */
+ if ((gparts[i].time_bin != time_bin_inhibited) &&
+ (gparts[i].time_bin != time_bin_not_created) &&
+ (gparts[i].type == swift_type_dark_matter_background)) {
+
+ if (with_stf) vr_data_written[count] = vr_data[i];
+
+ gparts_written[count] = gparts[i];
+ count++;
+ }
+ }
+
+ /* Check that everything is fine */
+ if (count != Ngparts_written)
+ error("Collected the wrong number of g-particles (%zu vs. %zu expected)",
+ count, Ngparts_written);
+}
+
/**
* @brief Verify the io parameter file
*
@@ -1895,11 +2031,7 @@ void io_collect_gparts_to_write(
* @param N_total The total number of each particle type.
*/
void io_check_output_fields(const struct swift_params* params,
- const long long N_total[3]) {
-
- /* Copy N_total to array with length == 6 */
- const long long nr_total[swift_type_count] = {N_total[0], N_total[1], 0,
- 0, N_total[2], 0};
+ const long long N_total[swift_type_count]) {
/* Loop over all particle types to check the fields */
for (int ptype = 0; ptype < swift_type_count; ptype++) {
@@ -1908,7 +2040,7 @@ void io_check_output_fields(const struct swift_params* params,
struct io_props list[100];
/* Don't do anything if no particle of this kind */
- if (nr_total[ptype] == 0) continue;
+ if (N_total[ptype] == 0) continue;
/* Gather particle fields from the particle structures */
switch (ptype) {
@@ -1932,6 +2064,12 @@ void io_check_output_fields(const struct swift_params* params,
num_fields += velociraptor_write_gparts(NULL, list + num_fields);
break;
+ case swift_type_dark_matter_background:
+ darkmatter_write_particles(NULL, list, &num_fields);
+ num_fields += fof_write_gparts(NULL, list + num_fields);
+ num_fields += velociraptor_write_gparts(NULL, list + num_fields);
+ break;
+
case swift_type_stars:
stars_write_particles(NULL, list, &num_fields, /*with_cosmology=*/1);
num_fields += chemistry_write_sparticles(NULL, list + num_fields);
@@ -2045,6 +2183,12 @@ void io_write_output_field_parameter(const char* filename) {
num_fields += velociraptor_write_gparts(NULL, list + num_fields);
break;
+ case swift_type_dark_matter_background:
+ darkmatter_write_particles(NULL, list, &num_fields);
+ num_fields += fof_write_gparts(NULL, list + num_fields);
+ num_fields += velociraptor_write_gparts(NULL, list + num_fields);
+ break;
+
case swift_type_stars:
stars_write_particles(NULL, list, &num_fields, /*with_cosmology=*/1);
num_fields += chemistry_write_sparticles(NULL, list + num_fields);
diff --git a/src/common_io.h b/src/common_io.h
index 03437d1a38a09e59e6abd4b9ca762f181ffc3731..aa4b0c552f8a416a4aa42ecce5e9ae0b74dbc731 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -125,8 +125,18 @@ void io_collect_gparts_to_write(const struct gpart* restrict gparts,
struct velociraptor_gpart_data* vr_data_written,
const size_t Ngparts,
const size_t Ngparts_written, int with_stf);
+void io_collect_gparts_background_to_write(
+ const struct gpart* restrict gparts,
+ const struct velociraptor_gpart_data* vr_data,
+ struct gpart* restrict gparts_written,
+ struct velociraptor_gpart_data* vr_data_written, const size_t Ngparts,
+ const size_t Ngparts_written, int with_stf);
void io_prepare_dm_gparts(struct threadpool* tp, struct gpart* const gparts,
size_t Ndm);
+void io_prepare_dm_background_gparts(struct threadpool* tp,
+ struct gpart* const gparts, size_t Ndm);
+size_t io_count_dm_background_gparts(const struct gpart* const gparts,
+ size_t Ndm);
void io_duplicate_hydro_gparts(struct threadpool* tp, struct part* const parts,
struct gpart* const gparts, size_t Ngas,
size_t Ndm);
diff --git a/src/debug.c b/src/debug.c
index f8509af805fed90cbdbdfd58849cb7b6352062e9..8e3c581163bb404b8013790ef32549503c93fa38 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -73,6 +73,8 @@
#include "./gravity/Default/gravity_debug.h"
#elif defined(POTENTIAL_GRAVITY)
#include "./gravity/Potential/gravity_debug.h"
+#elif defined(MULTI_SOFTENING_GRAVITY)
+#include "./gravity/MultiSoftening/gravity_debug.h"
#else
#error "Invalid choice of gravity variant"
#endif
diff --git a/src/engine.c b/src/engine.c
index c3ace9f9f7d237708f57f5ddfd575b11d3c62ed1..462f38ba4b95b7b6322964d91e19285eb618d48f 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -4201,7 +4201,7 @@ void engine_do_reconstruct_multipoles_mapper(void *map_data, int num_elements,
if (c != NULL && c->nodeID == e->nodeID) {
/* Construct the multipoles in this cell hierarchy */
- cell_make_multipoles(c, e->ti_current);
+ cell_make_multipoles(c, e->ti_current, e->gravity_properties);
}
}
}
@@ -4408,19 +4408,26 @@ void engine_makeproxies(struct engine *e) {
sqrt(min_dist_centres2) - 2. * delta_CoM;
const double min_dist_CoM2 = min_dist_CoM * min_dist_CoM;
+ /* We also assume that the softening is negligible compared
+ to the cell size */
+ const double epsilon_i = 0.;
+ const double epsilon_j = 0.;
+
/* Are we beyond the distance where the truncated forces are 0
* but not too far such that M2L can be used? */
if (periodic) {
if ((min_dist_CoM2 < max_mesh_dist2) &&
(!gravity_M2L_accept(r_max, r_max, theta_crit2,
- min_dist_CoM2)))
+ min_dist_CoM2, epsilon_i,
+ epsilon_j)))
proxy_type |= (int)proxy_cell_type_gravity;
} else {
if (!gravity_M2L_accept(r_max, r_max, theta_crit2,
- min_dist_CoM2))
+ min_dist_CoM2, epsilon_i,
+ epsilon_j))
proxy_type |= (int)proxy_cell_type_gravity;
}
}
@@ -4883,6 +4890,7 @@ void engine_unpin(void) {
* @param Ngparts total number of gravity particles in the simulation.
* @param Nstars total number of star particles in the simulation.
* @param Nblackholes total number of black holes in the simulation.
+ * @param Nbackground_gparts Total number of background DM particles.
* @param policy The queuing policy to use.
* @param verbose Is this #engine talkative ?
* @param reparttype What type of repartition algorithm are we using ?
@@ -4904,8 +4912,8 @@ void engine_unpin(void) {
*/
void engine_init(struct engine *e, struct space *s, struct swift_params *params,
long long Ngas, long long Ngparts, long long Nstars,
- long long Nblackholes, int policy, int verbose,
- struct repartition *reparttype,
+ long long Nblackholes, long long Nbackground_gparts,
+ int policy, int verbose, struct repartition *reparttype,
const struct unit_system *internal_units,
const struct phys_const *physical_constants,
struct cosmology *cosmo, struct hydro_props *hydro,
@@ -4930,6 +4938,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
e->total_nr_gparts = Ngparts;
e->total_nr_sparts = Nstars;
e->total_nr_bparts = Nblackholes;
+ e->total_nr_DM_background_gparts = Nbackground_gparts;
e->proxy_ind = NULL;
e->nr_proxies = 0;
e->reparttype = reparttype;
@@ -6155,11 +6164,17 @@ void engine_recompute_displacement_constraint(struct engine *e) {
#endif
/* Get the counts of each particle types */
+ const long long total_nr_baryons =
+ e->total_nr_parts + e->total_nr_sparts + e->total_nr_bparts;
const long long total_nr_dm_gparts =
- e->total_nr_gparts - e->total_nr_parts - e->total_nr_sparts;
+ e->total_nr_gparts - e->total_nr_DM_background_gparts - total_nr_baryons;
float count_parts[swift_type_count] = {
- (float)e->total_nr_parts, (float)total_nr_dm_gparts, 0.f, 0.f,
- (float)e->total_nr_sparts, (float)e->total_nr_bparts};
+ (float)e->total_nr_parts,
+ (float)total_nr_dm_gparts,
+ (float)e->total_nr_DM_background_gparts,
+ 0.f,
+ (float)e->total_nr_sparts,
+ (float)e->total_nr_bparts};
/* Count of particles for the two species */
const float N_dm = count_parts[1];
@@ -6559,7 +6574,8 @@ void engine_fof(struct engine *e, const int dump_results,
/* Compute number of DM particles */
const long long total_nr_baryons =
e->total_nr_parts + e->total_nr_sparts + e->total_nr_bparts;
- const long long total_nr_dmparts = e->total_nr_gparts - total_nr_baryons;
+ const long long total_nr_dmparts =
+ e->total_nr_gparts - e->total_nr_DM_background_gparts - total_nr_baryons;
/* Initialise FOF parameters and allocate FOF arrays. */
fof_allocate(e->s, total_nr_dmparts, e->fof_properties);
diff --git a/src/engine.h b/src/engine.h
index dc5c2c9ddd9a843d3e677e9b56c47f7451515b2f..3c74f89dfcade93bb67d21eb6607cdd4cb80155d 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -248,6 +248,7 @@ struct engine {
long long total_nr_gparts;
long long total_nr_sparts;
long long total_nr_bparts;
+ long long total_nr_DM_background_gparts;
/* Total numbers of cells (top-level and sub-cells) in the system. */
long long total_nr_cells;
@@ -493,8 +494,8 @@ void engine_dump_snapshot(struct engine *e);
void engine_init_output_lists(struct engine *e, struct swift_params *params);
void engine_init(struct engine *e, struct space *s, struct swift_params *params,
long long Ngas, long long Ngparts, long long Nstars,
- long long Nblackholes, int policy, int verbose,
- struct repartition *reparttype,
+ long long Nblackholes, long long Nbackground_gparts,
+ int policy, int verbose, struct repartition *reparttype,
const struct unit_system *internal_units,
const struct phys_const *physical_constants,
struct cosmology *cosmo, struct hydro_props *hydro,
diff --git a/src/fof.c b/src/fof.c
index c7c9323a3ed0e5c52cf09c2644c3754e6126a24f..af530af6aadfd8a8c1ea8529330156f84278a534 100644
--- a/src/fof.c
+++ b/src/fof.c
@@ -1956,14 +1956,14 @@ void fof_dump_group_data(const struct fof_props *props,
: -1;
#ifdef WITH_MPI
fprintf(file, " %8zu %12zu %12e %12e %18lld %18lld\n",
- gparts[group_offset - node_offset].fof_data.group_id,
+ (size_t)gparts[group_offset - node_offset].fof_data.group_id,
group_size[group_offset - node_offset], group_mass[i],
max_part_density[i], max_part_density_index[i], part_id);
#else
fprintf(file, " %8zu %12zu %12e %12e %18lld %18lld\n",
- gparts[group_offset].fof_data.group_id, group_size[group_offset],
- group_mass[i], max_part_density[i], max_part_density_index[i],
- part_id);
+ (size_t)gparts[group_offset].fof_data.group_id,
+ group_size[group_offset], group_mass[i], max_part_density[i],
+ max_part_density_index[i], part_id);
#endif
}
diff --git a/src/gravity.h b/src/gravity.h
index 57cad8eba5f0772f85affd031a450c3ecb4dde59..4e255d2cbaaae012ab1282c82c8f6f15410e6b4b 100644
--- a/src/gravity.h
+++ b/src/gravity.h
@@ -34,6 +34,9 @@
#elif defined(POTENTIAL_GRAVITY)
#include "./gravity/Potential/gravity.h"
#define GRAVITY_IMPLEMENTATION "With potential calculation"
+#elif defined(MULTI_SOFTENING_GRAVITY)
+#include "./gravity/MultiSoftening/gravity.h"
+#define GRAVITY_IMPLEMENTATION "With per-particle softening"
#else
#error "Invalid choice of gravity variant"
#endif
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 6d1270c952100f3a25202fcdb22be09f9acaa8d9..a8ace72d2b2bb8cb92d2d386ae3ed5402c4f76e0 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -42,13 +42,16 @@ __attribute__((always_inline)) INLINE static float gravity_get_mass(
/**
* @brief Returns the current co-moving softening of a particle
*
+ * Note that in this basic gravity scheme, all particles have
+ * the same softening length.
+ *
* @param gp The particle of interest
* @param grav_props The global gravity properties.
*/
__attribute__((always_inline)) INLINE static float gravity_get_softening(
const struct gpart* gp, const struct gravity_props* restrict grav_props) {
- return grav_props->epsilon_cur;
+ return grav_props->epsilon_DM_cur;
}
/**
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 6fce3ddd512018e9ea4be21111c75904c77cb925..4967b53d15d0c1888c77ef51265be31f571009cc 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -166,8 +166,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
/* Compute the derivatives of the potential */
struct potential_derivatives_M2P d;
- potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, 0.f,
- &d);
+ potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h,
+ /*periodic=*/0, /*r_s_inv=*/0.f, &d);
/* 0th order contributions */
*f_x = m->M_000 * d.D_100;
@@ -271,7 +271,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
/* Compute the derivatives of the potential */
struct potential_derivatives_M2P d;
- potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1,
+ potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, /*periodic=*/1,
r_s_inv, &d);
/* 0th order contributions */
diff --git a/src/gravity/MultiSoftening/gravity.h b/src/gravity/MultiSoftening/gravity.h
new file mode 100644
index 0000000000000000000000000000000000000000..08282a038dc67c2740e3d5e6ba95a5a5b8a190ff
--- /dev/null
+++ b/src/gravity/MultiSoftening/gravity.h
@@ -0,0 +1,246 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2016 Tom Theuns (tom.theuns@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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_MULTI_SOFTENING_GRAVITY_H
+#define SWIFT_MULIT_SOFTENING_GRAVITY_H
+
+#include <float.h>
+
+/* Local includes. */
+#include "cosmology.h"
+#include "gravity_properties.h"
+#include "kernel_gravity.h"
+#include "minmax.h"
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param gp The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float gravity_get_mass(
+ const struct gpart* restrict gp) {
+
+ return gp->mass;
+}
+
+/**
+ * @brief Returns the current co-moving softening of a particle
+ *
+ * @param gp The particle of interest
+ * @param grav_props The global gravity properties.
+ */
+__attribute__((always_inline)) INLINE static float gravity_get_softening(
+ const struct gpart* gp, const struct gravity_props* restrict grav_props) {
+
+ switch (gp->type) {
+ case swift_type_dark_matter:
+ return grav_props->epsilon_DM_cur;
+ case swift_type_stars:
+ return grav_props->epsilon_baryon_cur;
+ case swift_type_gas:
+ return grav_props->epsilon_baryon_cur;
+ case swift_type_black_hole:
+ return grav_props->epsilon_baryon_cur;
+ case swift_type_dark_matter_background:
+ return grav_props->epsilon_background_fac * cbrtf(gp->mass);
+ default:
+ return 0.f;
+ }
+}
+
+/**
+ * @brief Add a contribution to this particle's potential.
+ *
+ * @param gp The particle.
+ * @param pot The contribution to add.
+ */
+__attribute__((always_inline)) INLINE static void
+gravity_add_comoving_potential(struct gpart* restrict gp, float pot) {
+
+ gp->potential += pot;
+}
+
+/**
+ * @brief Returns the comoving potential of a particle.
+ *
+ * @param gp The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_get_comoving_potential(const struct gpart* restrict gp) {
+
+ return gp->potential;
+}
+
+/**
+ * @brief Returns the physical potential of a particle
+ *
+ * @param gp The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_get_physical_potential(const struct gpart* restrict gp,
+ const struct cosmology* cosmo) {
+
+ return gp->potential * cosmo->a_inv;
+}
+
+/**
+ * @brief Computes the gravity time-step of a given particle due to self-gravity
+ *
+ * We use Gadget-2's type 0 time-step criterion.
+ *
+ * @param gp Pointer to the g-particle data.
+ * @param a_hydro The accelerations coming from the hydro scheme (can be 0).
+ * @param grav_props Constants used in the gravity scheme.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_compute_timestep_self(const struct gpart* const gp,
+ const float a_hydro[3],
+ const struct gravity_props* restrict grav_props,
+ const struct cosmology* cosmo) {
+
+ /* Get physical acceleration (gravity contribution) */
+ float a_phys_x = gp->a_grav[0] * cosmo->a_factor_grav_accel;
+ float a_phys_y = gp->a_grav[1] * cosmo->a_factor_grav_accel;
+ float a_phys_z = gp->a_grav[2] * cosmo->a_factor_grav_accel;
+
+ /* Get physical acceleration (hydro contribution) */
+ a_phys_x += a_hydro[0] * cosmo->a_factor_hydro_accel;
+ a_phys_y += a_hydro[1] * cosmo->a_factor_hydro_accel;
+ a_phys_z += a_hydro[2] * cosmo->a_factor_hydro_accel;
+
+ const float ac2 =
+ a_phys_x * a_phys_x + a_phys_y * a_phys_y + a_phys_z * a_phys_z;
+
+ const float ac_inv = (ac2 > 0.f) ? 1.f / sqrtf(ac2) : FLT_MAX;
+
+ const float epsilon = gravity_get_softening(gp, grav_props);
+
+ const float dt = sqrtf(2. * kernel_gravity_softening_plummer_equivalent_inv *
+ cosmo->a * grav_props->eta * epsilon * ac_inv);
+
+ return dt;
+}
+
+/**
+ * @brief Prepares a g-particle for the gravity calculation
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the variaous tasks
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void gravity_init_gpart(
+ struct gpart* gp) {
+
+ /* Zero the acceleration */
+ gp->a_grav[0] = 0.f;
+ gp->a_grav[1] = 0.f;
+ gp->a_grav[2] = 0.f;
+ gp->potential = 0.f;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ gp->potential_PM = 0.f;
+ gp->a_grav_PM[0] = 0.f;
+ gp->a_grav_PM[1] = 0.f;
+ gp->a_grav_PM[2] = 0.f;
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gp->num_interacted = 0;
+ gp->initialised = 1;
+#endif
+}
+
+/**
+ * @brief Finishes the gravity calculation.
+ *
+ * Multiplies the forces and accelerations by the appropiate constants.
+ * Applies cosmological correction for periodic BCs.
+ *
+ * No need to apply the potential normalisation correction for periodic
+ * BCs here since we do not compute the potential.
+ *
+ * @param gp The particle to act upon
+ * @param const_G Newton's constant in internal units.
+ * @param potential_normalisation Term to be added to all the particles.
+ * @param periodic Are we using periodic BCs?
+ */
+__attribute__((always_inline)) INLINE static void gravity_end_force(
+ struct gpart* gp, float const_G, const float potential_normalisation,
+ const int periodic) {
+
+ /* Apply the periodic correction to the peculiar potential */
+ if (periodic) gp->potential += potential_normalisation;
+
+ /* Let's get physical... */
+ gp->a_grav[0] *= const_G;
+ gp->a_grav[1] *= const_G;
+ gp->a_grav[2] *= const_G;
+ gp->potential *= const_G;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ gp->potential_PM *= const_G;
+ gp->a_grav_PM[0] *= const_G;
+ gp->a_grav_PM[1] *= const_G;
+ gp->a_grav_PM[2] *= const_G;
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gp->initialised = 0; /* Ready for next step */
+#endif
+}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param gp The particle to act upon
+ * @param dt The time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void gravity_kick_extra(
+ struct gpart* gp, float dt) {}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param gp The particle.
+ */
+__attribute__((always_inline)) INLINE static void
+gravity_reset_predicted_values(struct gpart* gp) {}
+
+/**
+ * @brief Initialises the g-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param gp The particle to act upon
+ * @param grav_props The global properties of the gravity calculation.
+ */
+__attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
+ struct gpart* gp, const struct gravity_props* grav_props) {
+
+ gp->time_bin = 0;
+
+ gravity_init_gpart(gp);
+}
+
+#endif /* SWIFT_MULTI_SOFTENING_GRAVITY_H */
diff --git a/src/gravity/MultiSoftening/gravity_debug.h b/src/gravity/MultiSoftening/gravity_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..7deea93f385d77871ed4bb28b881d7e93ce59436
--- /dev/null
+++ b/src/gravity/MultiSoftening/gravity_debug.h
@@ -0,0 +1,37 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_MULTI_SOFTENING_GRAVITY_DEBUG_H
+#define SWIFT_MULTI_SOFTENING_GRAVITY_DEBUG_H
+
+__attribute__((always_inline)) INLINE static void gravity_debug_particle(
+ const struct gpart* p) {
+ printf(
+ "mass=%.3e time_bin=%d\n"
+ "x=[%.5e,%.5e,%.5e], v_full=[%.5e,%.5e,%.5e], type=%d, "
+ "a=[%.5e,%.5e,%.5e], pot=%.5e\n",
+ p->mass, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0],
+ p->v_full[1], p->v_full[2], (int)p->type, p->a_grav[0], p->a_grav[1],
+ p->a_grav[2], p->potential);
+#ifdef SWIFT_DEBUG_CHECKS
+ printf("num_interacted=%lld ti_drift=%lld ti_kick=%lld\n", p->num_interacted,
+ p->ti_drift, p->ti_kick);
+#endif
+}
+
+#endif /* SWIFT_MULTI_SOFTENING_GRAVITY_DEBUG_H */
diff --git a/src/gravity/MultiSoftening/gravity_iact.h b/src/gravity/MultiSoftening/gravity_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..8ab3a390010d2a2f72d0a056e399f7720e5dc464
--- /dev/null
+++ b/src/gravity/MultiSoftening/gravity_iact.h
@@ -0,0 +1,350 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_MULTI_SOFTENING_GRAVITY_IACT_H
+#define SWIFT_MULTI_SOFTENING_GRAVITY_IACT_H
+
+/* Includes. */
+#include "kernel_gravity.h"
+#include "kernel_long_gravity.h"
+#include "multipole.h"
+
+/* Standard headers */
+#include <float.h>
+
+/**
+ * @brief Computes the intensity of the force at a point generated by a
+ * point-mass.
+ *
+ * The returned quantity needs to be multiplied by the distance vector to obtain
+ * the force vector.
+ *
+ * @param r2 Square of the distance to the point-mass.
+ * @param h2 Square of the softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param h_inv3 Cube of the inverse of the softening length.
+ * @param mass Mass of the point-mass.
+ * @param f_ij (return) The force intensity.
+ * @param pot_ij (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
+ const float r2, const float h2, const float h_inv, const float h_inv3,
+ const float mass, float *f_ij, float *pot_ij) {
+
+ /* Get the inverse distance */
+ const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
+
+ /* Should we soften ? */
+ if (r2 >= h2) {
+
+ /* Get Newtonian gravity */
+ *f_ij = mass * r_inv * r_inv * r_inv;
+ *pot_ij = -mass * r_inv;
+
+ } else {
+
+ const float r = r2 * r_inv;
+ const float ui = r * h_inv;
+
+ float W_f_ij, W_pot_ij;
+ kernel_grav_force_eval(ui, &W_f_ij);
+ kernel_grav_pot_eval(ui, &W_pot_ij);
+
+ /* Get softened gravity */
+ *f_ij = mass * h_inv3 * W_f_ij;
+ *pot_ij = mass * h_inv * W_pot_ij;
+ }
+}
+
+/**
+ * @brief Computes the intensity of the force at a point generated by a
+ * point-mass truncated for long-distance periodicity.
+ *
+ * The returned quantity needs to be multiplied by the distance vector to obtain
+ * the force vector.
+ *
+ * @param r2 Square of the distance to the point-mass.
+ * @param h2 Square of the softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param h_inv3 Cube of the inverse of the softening length.
+ * @param mass Mass of the point-mass.
+ * @param r_s_inv Inverse of the mesh smoothing scale.
+ * @param f_ij (return) The force intensity.
+ * @param pot_ij (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
+ const float r2, const float h2, const float h_inv, const float h_inv3,
+ const float mass, const float r_s_inv, float *f_ij, float *pot_ij) {
+
+ /* Get the inverse distance */
+ const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
+ const float r = r2 * r_inv;
+
+ /* Should we soften ? */
+ if (r2 >= h2) {
+
+ /* Get Newtonian gravity */
+ *f_ij = mass * r_inv * r_inv * r_inv;
+ *pot_ij = -mass * r_inv;
+
+ } else {
+
+ const float ui = r * h_inv;
+ float W_f_ij, W_pot_ij;
+
+ kernel_grav_force_eval(ui, &W_f_ij);
+ kernel_grav_pot_eval(ui, &W_pot_ij);
+
+ /* Get softened gravity */
+ *f_ij = mass * h_inv3 * W_f_ij;
+ *pot_ij = mass * h_inv * W_pot_ij;
+ }
+
+ /* Get long-range correction */
+ const float u_lr = r * r_s_inv;
+ float corr_f_lr, corr_pot_lr;
+ kernel_long_grav_force_eval(u_lr, &corr_f_lr);
+ kernel_long_grav_pot_eval(u_lr, &corr_pot_lr);
+ *f_ij *= corr_f_lr;
+ *pot_ij *= corr_pot_lr;
+}
+
+/**
+ * @brief Computes the forces at a point generated by a multipole.
+ *
+ * This assumes M_100 == M_010 == M_001 == 0.
+ * This uses the quadrupole and trace of the octupole terms only and defaults to
+ * the monopole if the code is compiled with low-order gravity only.
+ *
+ * @param r_x x-component of the distance vector to the multipole.
+ * @param r_y y-component of the distance vector to the multipole.
+ * @param r_z z-component of the distance vector to the multipole.
+ * @param r2 Square of the distance vector to the multipole.
+ * @param h The softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param m The multipole.
+ * @param f_x (return) The x-component of the acceleration.
+ * @param f_y (return) The y-component of the acceleration.
+ * @param f_z (return) The z-component of the acceleration.
+ * @param pot (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
+ const float r_x, const float r_y, const float r_z, const float r2,
+ const float h, const float h_inv, const struct multipole *m, float *f_x,
+ float *f_y, float *f_z, float *pot) {
+
+/* In the case where the order is < 3, then there is only a monopole term left.
+ * We can default to the normal P-P interaction with the mass of the multipole
+ * and its CoM as the "particle" property */
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
+
+ float f_ij, pot_ij;
+ runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
+ &f_ij, &pot_ij);
+ *f_x = f_ij * r_x;
+ *f_y = f_ij * r_y;
+ *f_z = f_ij * r_z;
+ *pot = pot_ij;
+
+#else
+
+ /* Get the inverse distance */
+ const float r_inv = 1.f / sqrtf(r2);
+
+ /* Compute the derivatives of the potential */
+ struct potential_derivatives_M2P d;
+ potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h,
+ /*periodic=*/0, /*r_s_inv=*/0.f, &d);
+
+ /* 0th order contributions */
+ *f_x = m->M_000 * d.D_100;
+ *f_y = m->M_000 * d.D_010;
+ *f_z = m->M_000 * d.D_001;
+ *pot = m->M_000 * d.D_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
+ /* 1st order contributions */
+
+ /* 1st order contributions are all 0 since the dipole is 0 */
+
+ /* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */
+ /* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */
+ /* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
+ /* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
+
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* 2nd order contributions */
+ *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
+ m->M_101 * d.D_201 + m->M_110 * d.D_210 + m->M_200 * d.D_300;
+ *f_y += m->M_002 * d.D_012 + m->M_011 * d.D_021 + m->M_020 * d.D_030 +
+ m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210;
+ *f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 +
+ m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
+ *pot += m->M_002 * d.D_002 + m->M_011 * d.D_011 + m->M_020 * d.D_020 +
+ m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
+
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+
+ /* 3rd order contributions */
+ *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
+ m->M_030 * d.D_130 + m->M_102 * d.D_202 + m->M_111 * d.D_211 +
+ m->M_120 * d.D_220 + m->M_201 * d.D_301 + m->M_210 * d.D_310 +
+ m->M_300 * d.D_400;
+ *f_y += m->M_003 * d.D_013 + m->M_012 * d.D_022 + m->M_021 * d.D_031 +
+ m->M_030 * d.D_040 + m->M_102 * d.D_112 + m->M_111 * d.D_121 +
+ m->M_120 * d.D_130 + m->M_201 * d.D_211 + m->M_210 * d.D_220 +
+ m->M_300 * d.D_310;
+ *f_z += m->M_003 * d.D_004 + m->M_012 * d.D_013 + m->M_021 * d.D_022 +
+ m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 +
+ m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
+ m->M_300 * d.D_301;
+ *pot += m->M_003 * d.D_003 + m->M_012 * d.D_012 + m->M_021 * d.D_021 +
+ m->M_030 * d.D_030 + m->M_102 * d.D_102 + m->M_111 * d.D_111 +
+ m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
+ m->M_300 * d.D_300;
+
+#endif
+
+ /* Take care of the the sign convention */
+ *f_x *= -1.f;
+ *f_y *= -1.f;
+ *f_z *= -1.f;
+ *pot *= -1.f;
+#endif
+}
+
+/**
+ * @brief Computes the forces at a point generated by a multipole, truncated for
+ * long-range periodicity.
+ *
+ * This assumes M_100 == M_010 == M_001 == 0.
+ * This uses the quadrupole term and trace of the octupole terms only and
+ * defaults to the monopole if the code is compiled with low-order gravity only.
+ *
+ * @param r_x x-component of the distance vector to the multipole.
+ * @param r_y y-component of the distance vector to the multipole.
+ * @param r_z z-component of the distance vector to the multipole.
+ * @param r2 Square of the distance vector to the multipole.
+ * @param h The softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param r_s_inv The inverse of the gravity mesh-smoothing scale.
+ * @param m The multipole.
+ * @param f_x (return) The x-component of the acceleration.
+ * @param f_y (return) The y-component of the acceleration.
+ * @param f_z (return) The z-component of the acceleration.
+ * @param pot (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
+ const float r_x, const float r_y, const float r_z, const float r2,
+ const float h, const float h_inv, const float r_s_inv,
+ const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) {
+
+/* In the case where the order is < 3, then there is only a monopole term left.
+ * We can default to the normal P-P interaction with the mass of the multipole
+ * and its CoM as the "particle" property */
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
+
+ float f_ij, pot_ij;
+ runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
+ m->M_000, r_s_inv, &f_ij, &pot_ij);
+ *f_x = f_ij * r_x;
+ *f_y = f_ij * r_y;
+ *f_z = f_ij * r_z;
+ *pot = -pot_ij;
+
+#else
+
+ /* Get the inverse distance */
+ const float r_inv = 1.f / sqrtf(r2);
+
+ /* Compute the derivatives of the potential */
+ struct potential_derivatives_M2P d;
+ potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, /*periodic=*/1,
+ r_s_inv, &d);
+
+ /* 0th order contributions */
+ *f_x = m->M_000 * d.D_100;
+ *f_y = m->M_000 * d.D_010;
+ *f_z = m->M_000 * d.D_001;
+ *pot = m->M_000 * d.D_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
+ /* 1st order contributions */
+
+ /* 1st order contributions are all 0 since the dipole is 0 */
+
+ /* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */
+ /* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */
+ /* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
+ /* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
+
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* 2nd order contributions */
+ *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
+ m->M_101 * d.D_201 + m->M_110 * d.D_210 + m->M_200 * d.D_300;
+ *f_y += m->M_002 * d.D_012 + m->M_011 * d.D_021 + m->M_020 * d.D_030 +
+ m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210;
+ *f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 +
+ m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
+ *pot += m->M_002 * d.D_002 + m->M_011 * d.D_011 + m->M_020 * d.D_020 +
+ m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
+
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+
+ /* 3rd order contributions */
+ *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
+ m->M_030 * d.D_130 + m->M_102 * d.D_202 + m->M_111 * d.D_211 +
+ m->M_120 * d.D_220 + m->M_201 * d.D_301 + m->M_210 * d.D_310 +
+ m->M_300 * d.D_400;
+ *f_y += m->M_003 * d.D_013 + m->M_012 * d.D_022 + m->M_021 * d.D_031 +
+ m->M_030 * d.D_040 + m->M_102 * d.D_112 + m->M_111 * d.D_121 +
+ m->M_120 * d.D_130 + m->M_201 * d.D_211 + m->M_210 * d.D_220 +
+ m->M_300 * d.D_310;
+ *f_z += m->M_003 * d.D_004 + m->M_012 * d.D_013 + m->M_021 * d.D_022 +
+ m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 +
+ m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
+ m->M_300 * d.D_301;
+ *pot += m->M_003 * d.D_003 + m->M_012 * d.D_012 + m->M_021 * d.D_021 +
+ m->M_030 * d.D_030 + m->M_102 * d.D_102 + m->M_111 * d.D_111 +
+ m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
+ m->M_300 * d.D_300;
+
+#endif
+
+ /* Take care of the the sign convention */
+ *f_x *= -1.f;
+ *f_y *= -1.f;
+ *f_z *= -1.f;
+ *pot *= -1.f;
+#endif
+}
+
+#endif /* SWIFT_MULTI_SOFTENING_GRAVITY_IACT_H */
diff --git a/src/gravity/MultiSoftening/gravity_io.h b/src/gravity/MultiSoftening/gravity_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..412a270afef5d3e7ab75530184782479aaa46b82
--- /dev/null
+++ b/src/gravity/MultiSoftening/gravity_io.h
@@ -0,0 +1,138 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_MULTI_SOFTENING_GRAVITY_IO_H
+#define SWIFT_MULTI_SOFTENING_GRAVITY_IO_H
+
+#include "io_properties.h"
+
+INLINE static void convert_gpart_pos(const struct engine* e,
+ const struct gpart* gp, double* ret) {
+
+ if (e->s->periodic) {
+ ret[0] = box_wrap(gp->x[0], 0.0, e->s->dim[0]);
+ ret[1] = box_wrap(gp->x[1], 0.0, e->s->dim[1]);
+ ret[2] = box_wrap(gp->x[2], 0.0, e->s->dim[2]);
+ } else {
+ ret[0] = gp->x[0];
+ ret[1] = gp->x[1];
+ ret[2] = gp->x[2];
+ }
+}
+
+INLINE static void convert_gpart_vel(const struct engine* e,
+ const struct gpart* gp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, gp->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ ret[0] = gp->v_full[0] + gp->a_grav[0] * dt_kick_grav;
+ ret[1] = gp->v_full[1] + gp->a_grav[1] * dt_kick_grav;
+ ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav;
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a_inv;
+ ret[1] *= cosmo->a_inv;
+ ret[2] *= cosmo->a_inv;
+}
+
+INLINE static void convert_gpart_soft(const struct engine* e,
+ const struct gpart* gp, float* ret) {
+
+ ret[0] = kernel_gravity_softening_plummer_equivalent_inv *
+ gravity_get_softening(gp, e->gravity_properties);
+}
+
+/**
+ * @brief Specifies which g-particle fields to read from a dataset
+ *
+ * @param gparts The g-particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+INLINE static void darkmatter_read_particles(struct gpart* gparts,
+ struct io_props* list,
+ int* num_fields) {
+
+ /* Say how much we want to read */
+ *num_fields = 4;
+
+ /* List what we want to read */
+ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+ UNIT_CONV_LENGTH, gparts, x);
+ list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+ UNIT_CONV_SPEED, gparts, v_full);
+ list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+ gparts, mass);
+ list[3] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset);
+}
+
+/**
+ * @brief Specifies which g-particle fields to write to a dataset
+ *
+ * @param gparts The g-particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+INLINE static void darkmatter_write_particles(const struct gpart* gparts,
+ struct io_props* list,
+ int* num_fields) {
+
+ /* Say how much we want to write */
+ *num_fields = 5;
+
+ /* List what we want to write */
+ list[0] = io_make_output_field_convert_gpart(
+ "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, gparts,
+ convert_gpart_pos, "Co-moving position of the particles");
+
+ list[1] = io_make_output_field_convert_gpart(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, gparts, convert_gpart_vel,
+ "Peculiar velocities of the stars. This is a * dx/dt where x is the "
+ "co-moving position of the particles.");
+
+ list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f,
+ gparts, mass, "Masses of the particles");
+
+ list[3] = io_make_output_field(
+ "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, gparts,
+ id_or_neg_offset, "Unique ID of the particles");
+
+ list[4] = io_make_output_field_convert_gpart(
+ "Softenings", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, gparts, convert_gpart_soft,
+ "Co-moving Plummer-equivalent softening lengths of the particles.");
+}
+
+#endif /* SWIFT_MULTI_SOFTENING_GRAVITY_IO_H */
diff --git a/src/gravity/MultiSoftening/gravity_part.h b/src/gravity/MultiSoftening/gravity_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..de410fbbb988136af46c6f9085bdafee863b9f99
--- /dev/null
+++ b/src/gravity/MultiSoftening/gravity_part.h
@@ -0,0 +1,88 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_MULTI_SOFTENING_GRAVITY_PART_H
+#define SWIFT_MULTI_SOFTENING_GRAVITY_PART_H
+
+#include "fof_struct.h"
+
+/* Gravity particle. */
+struct gpart {
+
+ /*! Particle ID. If negative, it is the negative offset of the #part with
+ which this gpart is linked. */
+ long long id_or_neg_offset;
+
+ /*! Particle position. */
+ double x[3];
+
+ /*! Particle velocity. */
+ float v_full[3];
+
+ /*! Particle acceleration. */
+ float a_grav[3];
+
+ /*! Gravitational potential */
+ float potential;
+
+ /*! Particle mass. */
+ float mass;
+
+ /*! Particle FoF properties (group ID, group size, ...) */
+ struct fof_gpart_data fof_data;
+
+ /*! Time-step length */
+ timebin_t time_bin;
+
+ /*! Type of the #gpart (DM, gas, star, ...) */
+ enum part_type type;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Numer of gparts this gpart interacted with */
+ long long num_interacted;
+
+ /* Time of the last drift */
+ integertime_t ti_drift;
+
+ /* Time of the last kick */
+ integertime_t ti_kick;
+
+ /* Has this particle been initialised? */
+ int initialised;
+
+#endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+ /*! Acceleration taken from the mesh only */
+ float a_grav_PM[3];
+
+ /*! Potential taken from the mesh only */
+ float potential_PM;
+
+ /* Brute-force particle acceleration. */
+ double a_grav_exact[3];
+
+ /* Brute-force particle potential. */
+ double potential_exact;
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_MULTI_SOFTENING_GRAVITY_PART_H */
diff --git a/src/gravity/Potential/gravity.h b/src/gravity/Potential/gravity.h
index 7d38e9126f1a313b169092b080ebc312c4bbe1bc..3a3e9680016ef51c554691c6f43f8279cde8caf9 100644
--- a/src/gravity/Potential/gravity.h
+++ b/src/gravity/Potential/gravity.h
@@ -22,6 +22,7 @@
#include <float.h>
+/* Local includes. */
#include "cosmology.h"
#include "gravity_properties.h"
#include "kernel_gravity.h"
@@ -41,13 +42,16 @@ __attribute__((always_inline)) INLINE static float gravity_get_mass(
/**
* @brief Returns the current co-moving softening of a particle
*
+ * Note that in this basic gravity scheme, all particles have
+ * the same softening length.
+ *
* @param gp The particle of interest
* @param grav_props The global gravity properties.
*/
__attribute__((always_inline)) INLINE static float gravity_get_softening(
const struct gpart* gp, const struct gravity_props* restrict grav_props) {
- return grav_props->epsilon_cur;
+ return grav_props->epsilon_DM_cur;
}
/**
diff --git a/src/gravity/Potential/gravity_iact.h b/src/gravity/Potential/gravity_iact.h
index f2094f6ecd5b31b94ebfe7a64f42fbd289a0c81c..abd4b0718145b131b0c1ad79dfe05d89a55cbb28 100644
--- a/src/gravity/Potential/gravity_iact.h
+++ b/src/gravity/Potential/gravity_iact.h
@@ -169,8 +169,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
/* Compute the derivatives of the potential */
struct potential_derivatives_M2P d;
- potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, 0.f,
- &d);
+ potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, /*periodic=*/0,
+ /*r_s_inv=*/0.f, &d);
/* 0th order contributions */
*f_x = m->M_000 * d.D_100;
@@ -178,7 +178,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
*f_z = m->M_000 * d.D_001;
*pot = m->M_000 * d.D_000;
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 1st order contributions */
@@ -281,7 +281,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
/* Compute the derivatives of the potential */
struct potential_derivatives_M2P d;
- potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1,
+ potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, /*periodic=*/1,
r_s_inv, &d);
/* 0th order contributions */
@@ -290,7 +290,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
*f_z = m->M_000 * d.D_001;
*pot = m->M_000 * d.D_000;
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 1st order contributions */
diff --git a/src/gravity_cache.h b/src/gravity_cache.h
index 6d073db60171a9d30d6f397213849e4c3d5314a1..912961d2e06dc748039f54c72bd9b33fd2547df2 100644
--- a/src/gravity_cache.h
+++ b/src/gravity_cache.h
@@ -255,7 +255,8 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
const float r2 = dx * dx + dy * dy + dz * dz;
/* Check whether we can use the multipole instead of P-P */
- use_mpole[i] = allow_mpole && gravity_M2P_accept(r_max2, theta_crit2, r2);
+ use_mpole[i] =
+ allow_mpole && gravity_M2P_accept(r_max2, theta_crit2, r2, epsilon[i]);
}
#ifdef SWIFT_DEBUG_CHECKS
@@ -436,7 +437,7 @@ gravity_cache_populate_all_mpole(const timebin_t max_active_bin,
}
const float r2 = dx * dx + dy * dy + dz * dz;
- if (!gravity_M2P_accept(r_max2, theta_crit2, r2))
+ if (!gravity_M2P_accept(r_max2, theta_crit2, r2, epsilon[i]))
error("Using m-pole where the test fails");
#endif
}
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index 3dcffe1cc04c5e10d3b2353b7e21c532747c1475..2ec9bd58c9c7ef8703b6191d1346c3da74416101 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -194,7 +194,6 @@ potential_derivatives_flip_signs(struct potential_derivatives_M2L *pot) {
* @param r2 Square norm of distance vector
* @param r_inv Inverse norm of distance vector
* @param eps Softening length.
- * @param eps_inv Inverse of softening length.
* @param periodic Is the calculation periodic ?
* @param r_s_inv Inverse of the long-range gravity mesh smoothing length.
* @param pot (return) The structure containing all the derivatives.
@@ -203,10 +202,14 @@ __attribute__((always_inline)) INLINE static void
potential_derivatives_compute_M2L(const float r_x, const float r_y,
const float r_z, const float r2,
const float r_inv, const float eps,
- const float eps_inv, const int periodic,
- const float r_s_inv,
+ const int periodic, const float r_s_inv,
struct potential_derivatives_M2L *pot) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (r2 < 0.99f * eps * eps)
+ error("Computing M2L derivatives below softening length");
+#endif
+
float Dt_1;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
float Dt_3;
@@ -224,8 +227,8 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
float Dt_11;
#endif
- /* Un-softened un-truncated case (Newtonian potential) */
- if (!periodic && r2 > eps * eps) {
+ /* Un-truncated case (Newtonian potential) */
+ if (!periodic) {
Dt_1 = r_inv;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
@@ -248,8 +251,8 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
#error "Missing implementation for order >5"
#endif
- /* Un-softened truncated case */
- } else if (periodic && r2 > eps * eps) {
+ /* Truncated case */
+ } else {
/* Get the derivatives of the truncated potential */
const float r = r2 * r_inv;
@@ -291,38 +294,6 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
-#endif
-
- /* Softened case */
- } else {
- const float r = r2 * r_inv;
- const float u = r * eps_inv;
- const float u_inv = r_inv * eps;
-
- Dt_1 = eps_inv * D_soft_1(u, u_inv);
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
- const float eps_inv2 = eps_inv * eps_inv;
- const float eps_inv3 = eps_inv * eps_inv2;
- Dt_3 = -eps_inv3 * D_soft_3(u, u_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
- const float eps_inv5 = eps_inv3 * eps_inv2;
- Dt_5 = eps_inv5 * D_soft_5(u, u_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
- const float eps_inv7 = eps_inv5 * eps_inv2;
- Dt_7 = -eps_inv7 * D_soft_7(u, u_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
- const float eps_inv9 = eps_inv7 * eps_inv2;
- Dt_9 = eps_inv9 * D_soft_9(u, u_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
- const float eps_inv11 = eps_inv9 * eps_inv2;
- Dt_11 = -eps_inv11 * D_soft_11(u, u_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
-#error "Missing implementation for order >5"
#endif
}
@@ -450,7 +421,6 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
* @param r2 Square norm of distance vector
* @param r_inv Inverse norm of distance vector
* @param eps Softening length.
- * @param eps_inv Inverse of softening length.
* @param periodic Is the calculation using periodic BCs?
* @param r_s_inv The inverse of the gravity mesh-smoothing scale.
* @param pot (return) The structure containing all the derivatives.
@@ -459,10 +429,14 @@ __attribute__((always_inline)) INLINE static void
potential_derivatives_compute_M2P(const float r_x, const float r_y,
const float r_z, const float r2,
const float r_inv, const float eps,
- const float eps_inv, const int periodic,
- const float r_s_inv,
+ const int periodic, const float r_s_inv,
struct potential_derivatives_M2P *pot) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (r2 < 0.99f * eps * eps)
+ error("Computing M2P derivatives below softening length");
+#endif
+
float Dt_1;
float Dt_3;
float Dt_5;
@@ -471,8 +445,8 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
float Dt_9;
#endif
- /* Un-softened un-truncated case (Newtonian potential) */
- if (!periodic && r2 > eps * eps) {
+ /* Un-truncated case (Newtonian potential) */
+ if (!periodic) {
const float r_inv2 = r_inv * r_inv;
@@ -484,8 +458,8 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
Dt_9 = -7.f * Dt_7 * r_inv2; /* -105 / r^9 */
#endif
- /* Un-softened truncated case */
- } else if (periodic && r2 > eps * eps) {
+ /* Truncated case */
+ } else if (periodic) {
/* Get the derivatives of the truncated potential */
const float r = r2 * r_inv;
@@ -512,30 +486,6 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
45.f * r * r * d.chi_2 - 105.f * r * d.chi_1 + 105.f * d.chi_0) *
r_inv9;
#endif
-
- /* Softened case */
- } else {
-
- const float r = r2 * r_inv;
- const float u = r * eps_inv;
- const float u_inv = r_inv * eps;
- const float eps_inv2 = eps_inv * eps_inv;
-
- Dt_1 = eps_inv * D_soft_1(u, u_inv);
-
- const float eps_inv3 = eps_inv * eps_inv2;
- Dt_3 = -eps_inv3 * D_soft_3(u, u_inv);
-
- const float eps_inv5 = eps_inv3 * eps_inv2;
- Dt_5 = eps_inv5 * D_soft_5(u, u_inv);
-
- const float eps_inv7 = eps_inv5 * eps_inv2;
- Dt_7 = -eps_inv7 * D_soft_7(u, u_inv);
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
- const float eps_inv9 = eps_inv7 * eps_inv2;
- Dt_9 = eps_inv9 * D_soft_9(u, u_inv);
-#endif
}
/* Compute some powers of r_x, r_y and r_z */
diff --git a/src/gravity_iact.h b/src/gravity_iact.h
index 2aaf7219d0851058cba92d4686ed989572dbcaa4..4c8d707fe33d0955809f8cb394d75c3a2f580fb7 100644
--- a/src/gravity_iact.h
+++ b/src/gravity_iact.h
@@ -32,6 +32,8 @@
#include "./gravity/Default/gravity_iact.h"
#elif defined(POTENTIAL_GRAVITY)
#include "./gravity/Potential/gravity_iact.h"
+#elif defined(MULTI_SOFTENING_GRAVITY)
+#include "./gravity/MultiSoftening/gravity_iact.h"
#else
#error "Invalid choice of gravity variant"
#endif
diff --git a/src/gravity_io.h b/src/gravity_io.h
index 752ee906081d706865e62bae7c8f505e9ca64347..44f392ee06d43228f7aa46faa5c7d9015fde790b 100644
--- a/src/gravity_io.h
+++ b/src/gravity_io.h
@@ -30,6 +30,8 @@
#include "./gravity/Default/gravity_io.h"
#elif defined(POTENTIAL_GRAVITY)
#include "./gravity/Potential/gravity_io.h"
+#elif defined(MULTI_SOFTENING_GRAVITY)
+#include "./gravity/MultiSoftening/gravity_io.h"
#else
#error "Invalid choice of gravity variant"
#endif
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index 8d3dfe9c4683842b905959ab3c688cc1f23309e9..d1654a1e7830d31c20731f130389a04d69b148b3 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -40,8 +40,9 @@
void gravity_props_init(struct gravity_props *p, struct swift_params *params,
const struct phys_const *phys_const,
- const struct cosmology *cosmo, int with_cosmology,
- int periodic) {
+ const struct cosmology *cosmo, const int with_cosmology,
+ const int has_baryons, const int has_DM,
+ const int is_zoom_simulation, const int periodic) {
/* Tree updates */
p->rebuild_frequency =
@@ -89,14 +90,58 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
/* Softening parameters */
if (with_cosmology) {
- p->epsilon_comoving =
- parser_get_param_double(params, "Gravity:comoving_softening");
- p->epsilon_max_physical =
- parser_get_param_double(params, "Gravity:max_physical_softening");
+
+ if (has_DM) {
+ /* Maximal physical softening taken straight from the parameter file */
+ p->epsilon_DM_max_physical =
+ parser_get_param_double(params, "Gravity:max_physical_DM_softening");
+
+ /* Co-moving softenings taken straight from the parameter file */
+ p->epsilon_DM_comoving =
+ parser_get_param_double(params, "Gravity:comoving_DM_softening");
+ }
+
+ if (has_baryons) {
+ /* Maximal physical softening taken straight from the parameter file */
+ p->epsilon_baryon_max_physical = parser_get_param_double(
+ params, "Gravity:max_physical_baryon_softening");
+
+ /* Co-moving softenings taken straight from the parameter file */
+ p->epsilon_baryon_comoving =
+ parser_get_param_double(params, "Gravity:comoving_baryon_softening");
+ }
+
+ if (is_zoom_simulation) {
+
+ /* Compute the comoving softening length for background particles as
+ * a fraction of the mean inter-particle density of the background DM
+ * particles Since they have variable masses the mass factor will be
+ * multiplied in later on. Note that we already multiply in the conversion
+ * from Plummer -> real softening length */
+ const double ratio_background =
+ parser_get_param_double(params, "Gravity:softening_ratio_background");
+
+ const double mean_matter_density =
+ cosmo->Omega_m * cosmo->critical_density_0;
+
+ p->epsilon_background_fac = kernel_gravity_softening_plummer_equivalent *
+ ratio_background *
+ cbrt(1. / mean_matter_density);
+ }
+
} else {
- p->epsilon_max_physical =
- parser_get_param_double(params, "Gravity:max_physical_softening");
- p->epsilon_comoving = p->epsilon_max_physical;
+
+ if (has_DM) {
+ p->epsilon_DM_max_physical =
+ parser_get_param_double(params, "Gravity:max_physical_DM_softening");
+ }
+ if (has_baryons) {
+ p->epsilon_baryon_max_physical = parser_get_param_double(
+ params, "Gravity:max_physical_baryon_softening");
+ }
+
+ p->epsilon_DM_comoving = p->epsilon_DM_max_physical;
+ p->epsilon_baryon_comoving = p->epsilon_baryon_max_physical;
}
/* Copy over the gravitational constant */
@@ -109,23 +154,26 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
void gravity_props_update(struct gravity_props *p,
const struct cosmology *cosmo) {
- /* Current softening lengths */
- double softening;
- if (p->epsilon_comoving * cosmo->a > p->epsilon_max_physical)
- softening = p->epsilon_max_physical / cosmo->a;
+ /* Current softening length for the high-res. DM particles. */
+ double DM_softening, baryon_softening;
+ if (p->epsilon_DM_comoving * cosmo->a > p->epsilon_DM_max_physical)
+ DM_softening = p->epsilon_DM_max_physical / cosmo->a;
else
- softening = p->epsilon_comoving;
+ DM_softening = p->epsilon_DM_comoving;
+
+ /* Current softening length for the high-res. baryon particles. */
+ if (p->epsilon_baryon_comoving * cosmo->a > p->epsilon_baryon_max_physical)
+ baryon_softening = p->epsilon_baryon_max_physical / cosmo->a;
+ else
+ baryon_softening = p->epsilon_baryon_comoving;
/* Plummer equivalent -> internal */
- softening *= kernel_gravity_softening_plummer_equivalent;
+ DM_softening *= kernel_gravity_softening_plummer_equivalent;
+ baryon_softening *= kernel_gravity_softening_plummer_equivalent;
/* Store things */
- p->epsilon_cur = softening;
-
- /* Other factors */
- p->epsilon_cur2 = softening * softening;
- p->epsilon_cur_inv = 1. / softening;
- p->epsilon_cur_inv3 = 1. / (softening * softening * softening);
+ p->epsilon_DM_cur = DM_softening;
+ p->epsilon_baryon_cur = baryon_softening;
}
void gravity_props_print(const struct gravity_props *p) {
@@ -143,16 +191,31 @@ void gravity_props_print(const struct gravity_props *p) {
kernel_gravity_softening_name);
message(
- "Self-gravity comoving softening: epsilon=%.4f (Plummer equivalent: "
- "%.4f)",
- p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent,
- p->epsilon_comoving);
+ "Self-gravity DM comoving softening: epsilon=%.6f (Plummer equivalent: "
+ "%.6f)",
+ p->epsilon_DM_comoving * kernel_gravity_softening_plummer_equivalent,
+ p->epsilon_DM_comoving);
message(
- "Self-gravity maximal physical softening: epsilon=%.4f (Plummer "
- "equivalent: %.4f)",
- p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent,
- p->epsilon_max_physical);
+ "Self-gravity DM maximal physical softening: epsilon=%.6f (Plummer "
+ "equivalent: %.6f)",
+ p->epsilon_DM_max_physical * kernel_gravity_softening_plummer_equivalent,
+ p->epsilon_DM_max_physical);
+
+ message(
+ "Self-gravity baryon comoving softening: epsilon=%.6f (Plummer "
+ "equivalent: "
+ "%.6f)",
+ p->epsilon_baryon_comoving * kernel_gravity_softening_plummer_equivalent,
+ p->epsilon_baryon_comoving);
+
+ message(
+ "Self-gravity baryon maximal physical softening: epsilon=%.6f "
+ "(Plummer "
+ "equivalent: %.6f)",
+ p->epsilon_baryon_max_physical *
+ kernel_gravity_softening_plummer_equivalent,
+ p->epsilon_baryon_max_physical);
message("Self-gravity mesh side-length: N=%d", p->mesh_size);
message("Self-gravity mesh smoothing-scale: a_smooth=%f", p->a_smooth);
@@ -174,20 +237,40 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
io_write_attribute_s(h_grpgrav, "Softening style",
kernel_gravity_softening_name);
+
io_write_attribute_f(
- h_grpgrav, "Comoving softening length [internal units]",
- p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent);
+ h_grpgrav, "Comoving DM softening length [internal units]",
+ p->epsilon_DM_comoving * kernel_gravity_softening_plummer_equivalent);
io_write_attribute_f(
h_grpgrav,
- "Comoving Softening length (Plummer equivalent) [internal units]",
- p->epsilon_comoving);
+ "Comoving DM softening length (Plummer equivalent) [internal units]",
+ p->epsilon_DM_comoving);
+
io_write_attribute_f(
- h_grpgrav, "Maximal physical softening length [internal units]",
- p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent);
+ h_grpgrav, "Maximal physical DM softening length [internal units]",
+ p->epsilon_DM_max_physical * kernel_gravity_softening_plummer_equivalent);
io_write_attribute_f(h_grpgrav,
- "Maximal physical softening length (Plummer equivalent) "
- " [internal units]",
- p->epsilon_max_physical);
+ "Maximal physical DM softening length (Plummer "
+ "equivalent) [internal units]",
+ p->epsilon_DM_max_physical);
+
+ io_write_attribute_f(
+ h_grpgrav, "Comoving baryon softening length [internal units]",
+ p->epsilon_baryon_comoving * kernel_gravity_softening_plummer_equivalent);
+ io_write_attribute_f(
+ h_grpgrav,
+ "Comoving baryon softening length (Plummer equivalent) [internal units]",
+ p->epsilon_baryon_comoving);
+
+ io_write_attribute_f(
+ h_grpgrav, "Maximal physical baryon softening length [internal units]",
+ p->epsilon_baryon_max_physical *
+ kernel_gravity_softening_plummer_equivalent);
+ io_write_attribute_f(h_grpgrav,
+ "Maximal physical baryon softening length (Plummer "
+ "equivalent) [internal units]",
+ p->epsilon_baryon_max_physical);
+
io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
io_write_attribute_s(h_grpgrav, "Scheme", GRAVITY_IMPLEMENTATION);
io_write_attribute_i(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
diff --git a/src/gravity_properties.h b/src/gravity_properties.h
index 79d72765e767351e29ffa30b615c1f7f4f98022c..d0a68ceb14d13f618dc628fb7e0a0dbf37d67cfb 100644
--- a/src/gravity_properties.h
+++ b/src/gravity_properties.h
@@ -40,25 +40,23 @@ struct swift_params;
*/
struct gravity_props {
- /*! Frequency of tree-rebuild in units of #gpart updates. */
- float rebuild_frequency;
+ /* -------------- Softening for the regular particles ---------------- */
- /*! Periodic long-range mesh side-length */
- int mesh_size;
+ /*! Co-moving softening length for high-res. DM particles at the current
+ * redshift. */
+ float epsilon_DM_cur;
- /*! Mesh smoothing scale in units of top-level cell size */
- float a_smooth;
+ /*! Co-moving softening length for high-res. baryon particles at the current
+ * redshift. */
+ float epsilon_baryon_cur;
- /*! Distance below which the truncated mesh force is Newtonian in units of
- * a_smooth */
- float r_cut_min_ratio;
+ /* -------------- Softening for the background DM -------------------- */
- /*! Distance above which the truncated mesh force is negligible in units of
- * a_smooth */
- float r_cut_max_ratio;
+ /*! Conversion factor from cbrt of particle mass to softening assuming
+ * a constant fraction of the mean inter-particle separation at that mass. */
+ float epsilon_background_fac;
- /*! Time integration dimensionless multiplier */
- float eta;
+ /* -------------- Properties of the FFM gravity ---------------------- */
/*! Tree opening angle (Multipole acceptance criterion) */
double theta_crit;
@@ -69,23 +67,49 @@ struct gravity_props {
/*! Inverse of opening angle */
double theta_crit_inv;
- /*! Comoving softening */
- double epsilon_comoving;
+ /* ------------- Properties of the softened gravity ------------------ */
+
+ /*! Co-moving softening length for for high-res. DM particles */
+ float epsilon_DM_comoving;
+
+ /*! Maximal softening length in physical coordinates for the high-res.
+ * DM particles */
+ float epsilon_DM_max_physical;
+
+ /*! Co-moving softening length for for high-res. baryon particles */
+ float epsilon_baryon_comoving;
+
+ /*! Maximal softening length in physical coordinates for the high-res.
+ * baryon particles */
+ float epsilon_baryon_max_physical;
- /*! Maxium physical softening */
- double epsilon_max_physical;
+ /*! Fraction of the mean inter particle separation corresponding to the
+ * co-moving softening length of the low-res. particles (DM + baryons) */
+ float mean_inter_particle_fraction_high_res;
- /*! Current softening length */
- float epsilon_cur;
+ /* ------------- Properties of the time integration ----------------- */
- /*! Square of current softening length */
- float epsilon_cur2;
+ /*! Frequency of tree-rebuild in units of #gpart updates. */
+ float rebuild_frequency;
+
+ /*! Time integration dimensionless multiplier */
+ float eta;
+
+ /* ------------- Properties of the mesh-based gravity ---------------- */
+
+ /*! Periodic long-range mesh side-length */
+ int mesh_size;
- /*! Inverse of current softening length */
- float epsilon_cur_inv;
+ /*! Mesh smoothing scale in units of top-level cell size */
+ float a_smooth;
- /*! Cube of the inverse of current oftening length */
- float epsilon_cur_inv3;
+ /*! Distance below which the truncated mesh force is Newtonian in units of
+ * a_smooth */
+ float r_cut_min_ratio;
+
+ /*! Distance above which the truncated mesh force is negligible in units of
+ * a_smooth */
+ float r_cut_max_ratio;
/*! Gravitational constant (in internal units, copied from the physical
* constants) */
@@ -95,8 +119,9 @@ struct gravity_props {
void gravity_props_print(const struct gravity_props *p);
void gravity_props_init(struct gravity_props *p, struct swift_params *params,
const struct phys_const *phys_const,
- const struct cosmology *cosmo, int with_cosmology,
- int periodic);
+ const struct cosmology *cosmo, const int with_cosmology,
+ const int has_baryons, const int has_DM,
+ const int is_zoom_simulation, const int periodic);
void gravity_props_update(struct gravity_props *p,
const struct cosmology *cosmo);
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index d272fb7381664d91f763119ae68043ce59a77e1f..d28cc1dff363517e72ebb3607ec5aa8121f828b2 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -343,7 +343,7 @@ void hydro_props_update(struct hydro_props *p, const struct gravity_props *gp,
* is a fixed fraction of the radius at which the softened forces
* recover a Newtonian behaviour (i.e. 2.8 * Plummer equivalent softening
* in the case of a cubic spline kernel). */
- p->h_min = p->h_min_ratio * gp->epsilon_cur / kernel_gamma;
+ p->h_min = p->h_min_ratio * gp->epsilon_baryon_cur / kernel_gamma;
}
/**
diff --git a/src/multipole.h b/src/multipole.h
index e867dfd4e2cc5c9fcd06d7d95dcf76a97689c2b3..bbe8c49db544425f4f35d57884c6036d3f83c905 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -117,6 +117,9 @@ struct multipole {
/*! Minimal velocity along each axis of all #gpart */
float min_delta_vel[3];
+ /*! Maximal co-moving softening of all the #gpart in the mulipole */
+ float max_softening;
+
/* 0th order term */
float M_000;
@@ -459,6 +462,7 @@ INLINE static void gravity_multipole_init(struct multipole *m) {
*/
INLINE static void gravity_multipole_print(const struct multipole *m) {
+ printf("eps_max = %12.5e\n", m->max_softening);
printf("Vel= [%12.5e %12.5e %12.5e]\n", m->vel[0], m->vel[1], m->vel[2]);
printf("-------------------------\n");
printf("M_000= %12.5e\n", m->M_000);
@@ -512,6 +516,9 @@ INLINE static void gravity_multipole_print(const struct multipole *m) {
INLINE static void gravity_multipole_add(struct multipole *restrict ma,
const struct multipole *restrict mb) {
+ /* Maximum of both softenings */
+ ma->max_softening = max(ma->max_softening, mb->max_softening);
+
/* Add 0th order term */
ma->M_000 += mb->M_000;
@@ -630,6 +637,14 @@ INLINE static int gravity_multipole_equal(const struct gravity_tensors *ga,
const double v2 = ma->vel[0] * ma->vel[0] + ma->vel[1] * ma->vel[1] +
ma->vel[2] * ma->vel[2];
+ /* Check maximal softening */
+ if (fabsf(ma->max_softening - mb->max_softening) /
+ fabsf(ma->max_softening + mb->max_softening) >
+ tolerance) {
+ message("max softening different!");
+ return 0;
+ }
+
/* Check bulk velocity (if non-zero and component > 1% of norm)*/
if (fabsf(ma->vel[0] + mb->vel[0]) > 1e-10 &&
(ma->vel[0] * ma->vel[0]) > 0.0001 * v2 &&
@@ -1022,11 +1037,14 @@ INLINE static int gravity_multipole_equal(const struct gravity_tensors *ga,
* @param multi The #multipole (content will be overwritten).
* @param gparts The #gpart.
* @param gcount The number of particles.
+ * @param grav_props The properties of the gravity scheme.
*/
INLINE static void gravity_P2M(struct gravity_tensors *multi,
- const struct gpart *gparts, int gcount) {
+ const struct gpart *gparts, const int gcount,
+ const struct gravity_props *const grav_props) {
/* Temporary variables */
+ float epsilon_max = 0.f;
double mass = 0.0;
double com[3] = {0.0, 0.0, 0.0};
double vel[3] = {0.f, 0.f, 0.f};
@@ -1034,12 +1052,14 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
/* Collect the particle data for CoM. */
for (int k = 0; k < gcount; k++) {
const double m = gparts[k].mass;
+ const float epsilon = gravity_get_softening(&gparts[k], grav_props);
#ifdef SWIFT_DEBUG_CHECKS
if (gparts[k].time_bin == time_bin_inhibited)
error("Inhibited particle in P2M. Should have been removed earlier.");
#endif
+ epsilon_max = max(epsilon_max, epsilon);
mass += m;
com[0] += gparts[k].x[0] * m;
com[1] += gparts[k].x[1] * m;
@@ -1203,6 +1223,7 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
#endif
/* Store the data on the multipole. */
+ multi->m_pole.max_softening = epsilon_max;
multi->m_pole.M_000 = mass;
multi->r_max = sqrt(r_max2);
multi->CoM[0] = com[0];
@@ -1316,6 +1337,9 @@ INLINE static void gravity_M2M(struct multipole *restrict m_a,
const struct multipole *restrict m_b,
const double pos_a[3], const double pos_b[3]) {
+ /* "shift" the softening */
+ m_a->max_softening = m_b->max_softening;
+
/* Shift 0th order term */
m_a->M_000 = m_b->M_000;
@@ -1964,8 +1988,7 @@ INLINE static void gravity_M2L_nonsym(
const int periodic, const double dim[3], const float rs_inv) {
/* Recover some constants */
- const float eps = props->epsilon_cur;
- const float eps_inv = props->epsilon_cur_inv;
+ const float eps = m_a->max_softening;
/* Compute distance vector */
float dx = (float)(pos_b[0] - pos_a[0]);
@@ -1985,8 +2008,8 @@ INLINE static void gravity_M2L_nonsym(
/* Compute all derivatives */
struct potential_derivatives_M2L pot;
- potential_derivatives_compute_M2L(dx, dy, dz, r2, r_inv, eps, eps_inv,
- periodic, rs_inv, &pot);
+ potential_derivatives_compute_M2L(dx, dy, dz, r2, r_inv, eps, periodic,
+ rs_inv, &pot);
/* Do the M2L tensor multiplication */
gravity_M2L_apply(l_b, m_a, &pot);
@@ -2015,8 +2038,7 @@ INLINE static void gravity_M2L_symmetric(
const float rs_inv) {
/* Recover some constants */
- const float eps = props->epsilon_cur;
- const float eps_inv = props->epsilon_cur_inv;
+ const float eps = max(m_a->max_softening, m_b->max_softening);
/* Compute distance vector */
float dx = (float)(pos_b[0] - pos_a[0]);
@@ -2036,8 +2058,8 @@ INLINE static void gravity_M2L_symmetric(
/* Compute all derivatives */
struct potential_derivatives_M2L pot;
- potential_derivatives_compute_M2L(dx, dy, dz, r2, r_inv, eps, eps_inv,
- periodic, rs_inv, &pot);
+ potential_derivatives_compute_M2L(dx, dy, dz, r2, r_inv, eps, periodic,
+ rs_inv, &pot);
/* Do the first M2L tensor multiplication */
gravity_M2L_apply(l_b, m_a, &pot);
@@ -2551,23 +2573,31 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
* We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
* Issue 1, pp.27-42, equation 10.
*
+ * We also additionally check that the distance between the multipoles
+ * is larger than the softening lengths (here the distance at which
+ * the gravity becomes Newtonian again, not the Plummer-equivalent quantity).
+ *
* @param r_crit_a The size of the multipole A.
* @param r_crit_b The size of the multipole B.
* @param theta_crit2 The square of the critical opening angle.
* @param r2 Square of the distance (periodically wrapped) between the
* multipoles.
+ * @param epsilon_a The maximal softening length of any particle in A.
+ * @param epsilon_b The maximal softening length of any particle in B.
*/
__attribute__((always_inline, const)) INLINE static int gravity_M2L_accept(
const double r_crit_a, const double r_crit_b, const double theta_crit2,
- const double r2) {
+ const double r2, const double epsilon_a, const double epsilon_b) {
const double size = r_crit_a + r_crit_b;
const double size2 = size * size;
+ const double epsilon_a2 = epsilon_a * epsilon_a;
+ const double epsilon_b2 = epsilon_b * epsilon_b;
// MATTHIEU: Make this mass-dependent ?
/* Multipole acceptance criterion (Dehnen 2002, eq.10) */
- return (r2 * theta_crit2 > size2);
+ return (r2 * theta_crit2 > size2) && (r2 > epsilon_a2) && (r2 > epsilon_b2);
}
/**
@@ -2577,18 +2607,24 @@ __attribute__((always_inline, const)) INLINE static int gravity_M2L_accept(
* We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
* Issue 1, pp.27-42, equation 10.
*
+ * We also additionally check that the distance between the particle and the
+ * multipole is larger than the softening length (here the distance at which
+ * the gravity becomes Newtonian again, not the Plummer-equivalent quantity).
+ *
* @param r_max2 The square of the size of the multipole.
* @param theta_crit2 The square of the critical opening angle.
* @param r2 Square of the distance (periodically wrapped) between the
* particle and the multipole.
+ * @param epsilon The softening length of the particle.
*/
__attribute__((always_inline, const)) INLINE static int gravity_M2P_accept(
- const float r_max2, const float theta_crit2, const float r2) {
+ const float r_max2, const float theta_crit2, const float r2,
+ const float epsilon) {
// MATTHIEU: Make this mass-dependent ?
/* Multipole acceptance criterion (Dehnen 2002, eq.10) */
- return (r2 * theta_crit2 > r_max2);
+ return (r2 * theta_crit2 > r_max2) && (r2 > epsilon * epsilon);
}
#endif /* SWIFT_MULTIPOLE_H */
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 929e0bb43ac15526bffe66a9563e38e5ea7af9a4..8ce21085e0dc0d1d126ed510262e14842ddf4019 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -704,12 +704,13 @@ void writeArray(struct engine* e, hid_t grp, char* fileName,
void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, struct bpart** bparts,
- size_t* Ngas, size_t* Ngparts, size_t* Nstars,
- size_t* Nblackholes, int* flag_entropy, int with_hydro,
- int with_gravity, int with_stars, int with_black_holes,
- int cleanup_h, int cleanup_sqrt_a, double h, double a,
- int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
- int n_threads, int dry_run) {
+ size_t* Ngas, size_t* Ngparts, size_t* Ngparts_background,
+ size_t* Nstars, size_t* Nblackholes, int* flag_entropy,
+ int with_hydro, int with_gravity, int with_stars,
+ int with_black_holes, int cleanup_h, int cleanup_sqrt_a,
+ double h, double a, int mpi_rank, int mpi_size,
+ MPI_Comm comm, MPI_Info info, int n_threads,
+ int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
@@ -721,9 +722,11 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
long long offset[swift_type_count] = {0};
int dimension = 3; /* Assume 3D if nothing is specified */
size_t Ndm = 0;
+ size_t Ndm_background = 0;
/* Initialise counters */
- *Ngas = 0, *Ngparts = 0, *Nstars = 0, *Nblackholes = 0;
+ *Ngas = 0, *Ngparts = 0, *Ngparts_background = 0, *Nstars = 0,
+ *Nblackholes = 0;
/* Open file */
/* message("Opening file '%s' as IC.", fileName); */
@@ -875,11 +878,14 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
/* Allocate memory to store gravity particles */
if (with_gravity) {
- Ndm = N[1];
+ Ndm = N[swift_type_dark_matter];
+ Ndm_background = N[swift_type_dark_matter_background];
*Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
N[swift_type_dark_matter] +
+ N[swift_type_dark_matter_background] +
(with_stars ? N[swift_type_stars] : 0) +
(with_black_holes ? N[swift_type_black_hole] : 0);
+ *Ngparts_background = Ndm_background;
if (swift_memalign("gparts", (void**)gparts, gpart_align,
*Ngparts * sizeof(struct gpart)) != 0)
error("Error while allocating memory for gravity particles");
@@ -928,6 +934,13 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
}
break;
+ case swift_type_dark_matter_background:
+ if (with_gravity) {
+ Nparticles = Ndm_background;
+ darkmatter_read_particles(*gparts + Ndm, list, &num_fields);
+ }
+ break;
+
case swift_type_stars:
if (with_stars) {
Nparticles = *Nstars;
@@ -968,17 +981,23 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
/* Prepare the DM particles */
io_prepare_dm_gparts(&tp, *gparts, Ndm);
+ /* Prepare the DM background particles */
+ io_prepare_dm_background_gparts(&tp, *gparts + Ndm, Ndm_background);
+
/* Duplicate the hydro particles into gparts */
- if (with_hydro) io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas, Ndm);
+ if (with_hydro)
+ io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas,
+ Ndm + Ndm_background);
/* Duplicate the stars particles into gparts */
if (with_stars)
- io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
+ io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars,
+ Ndm + Ndm_background + *Ngas);
/* Duplicate the stars particles into gparts */
if (with_black_holes)
io_duplicate_black_holes_gparts(&tp, *bparts, *gparts, *Nblackholes,
- Ndm + *Ngas + *Nstars);
+ Ndm + Ndm_background + *Ngas + *Nstars);
threadpool_clean(&tp);
}
@@ -1236,6 +1255,17 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
}
break;
+ case swift_type_dark_matter_background:
+ darkmatter_write_particles(gparts, list, &num_fields);
+ if (with_fof) {
+ num_fields += fof_write_gparts(gparts, list + num_fields);
+ }
+ if (with_stf) {
+ num_fields += velociraptor_write_gparts(e->s->gpart_group_data,
+ list + num_fields);
+ }
+ break;
+
case swift_type_stars:
stars_write_particles(sparts, list, &num_fields, with_cosmology);
num_fields += chemistry_write_sparticles(sparts, list + num_fields);
@@ -1329,6 +1359,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
const int with_cooling = e->policy & engine_policy_cooling;
const int with_temperature = e->policy & engine_policy_temperature;
const int with_fof = e->policy & engine_policy_fof;
+ const int with_DM_background = e->s->with_DM_background;
#ifdef HAVE_VELOCIRAPTOR
const int with_stf = (e->policy & engine_policy_structure_finding) &&
(e->s->gpart_group_data != NULL);
@@ -1344,6 +1375,11 @@ void write_output_parallel(struct engine* e, const char* baseName,
// const size_t Nbaryons = Ngas + Nstars;
// const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
+ size_t Ndm_background = 0;
+ if (with_DM_background) {
+ Ndm_background = io_count_dm_background_gparts(gparts, Ntot);
+ }
+
/* Number of particles that we will write */
const size_t Ntot_written =
e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
@@ -1356,10 +1392,11 @@ void write_output_parallel(struct engine* e, const char* baseName,
const size_t Nbaryons_written =
Ngas_written + Nstars_written + Nblackholes_written;
const size_t Ndm_written =
- Ntot_written > 0 ? Ntot_written - Nbaryons_written : 0;
+ Ntot_written > 0 ? Ntot_written - Nbaryons_written - Ndm_background : 0;
/* Compute offset in the file and total number of particles */
- size_t N[swift_type_count] = {Ngas_written, Ndm_written, 0, 0,
+ size_t N[swift_type_count] = {Ngas_written, Ndm_written,
+ Ndm_background, 0,
Nstars_written, Nblackholes_written};
long long N_total[swift_type_count] = {0};
long long offset[swift_type_count] = {0};
@@ -1665,6 +1702,43 @@ void write_output_parallel(struct engine* e, const char* baseName,
}
} break;
+ case swift_type_dark_matter_background: {
+
+ /* Ok, we need to fish out the particles we want */
+ Nparticles = Ndm_background;
+
+ /* Allocate temporary array */
+ if (swift_memalign("gparts_written", (void**)&gparts_written,
+ gpart_align,
+ Ndm_background * sizeof(struct gpart)) != 0)
+ error("Error while allocating temporart memory for gparts");
+
+ if (with_stf) {
+ if (swift_memalign(
+ "gpart_group_written", (void**)&gpart_group_data_written,
+ gpart_align,
+ Ndm_background * sizeof(struct velociraptor_gpart_data)) != 0)
+ error(
+ "Error while allocating temporart memory for gparts STF "
+ "data");
+ }
+
+ /* Collect the non-inhibited DM particles from gpart */
+ io_collect_gparts_background_to_write(
+ gparts, e->s->gpart_group_data, gparts_written,
+ gpart_group_data_written, Ntot, Ndm_background, with_stf);
+
+ /* Select the fields to write */
+ darkmatter_write_particles(gparts_written, list, &num_fields);
+ if (with_stf) {
+#ifdef HAVE_VELOCIRAPTOR
+ num_fields += velociraptor_write_gparts(gpart_group_data_written,
+ list + num_fields);
+#endif
+ }
+
+ } break;
+
case swift_type_stars: {
if (Nstars == Nstars_written) {
diff --git a/src/parallel_io.h b/src/parallel_io.h
index 6e6b52b74ebf166068f8dff110ea0ca0fd8464a2..465835013847347a907b3f81fd3611992070d0d7 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -36,12 +36,13 @@
void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, struct bpart** bparts,
- size_t* Ngas, size_t* Ngparts, size_t* Nsparts,
- size_t* Nbparts, int* flag_entropy, int with_hydro,
- int with_gravity, int with_stars, int with_black_holes,
- int cleanup_h, int cleanup_sqrt_a, double h, double a,
- int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
- int nr_threads, int dry_run);
+ size_t* Ngas, size_t* Ngparts, size_t* Ngparts_background,
+ size_t* Nsparts, size_t* Nbparts, int* flag_entropy,
+ int with_hydro, int with_gravity, int with_stars,
+ int with_black_holes, int cleanup_h, int cleanup_sqrt_a,
+ double h, double a, int mpi_rank, int mpi_size,
+ MPI_Comm comm, MPI_Info info, int nr_threads,
+ int dry_run);
void write_output_parallel(struct engine* e, const char* baseName,
const struct unit_system* internal_units,
diff --git a/src/part.c b/src/part.c
index 7d967caad18caf0c45d2a422fb4766d34f559d00..d8f57a8b3db3d5f29667dff162fb03d8b117d3a3 100644
--- a/src/part.c
+++ b/src/part.c
@@ -188,6 +188,15 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
error("DM gpart particle linked to something !");
}
+ /* We have a background DM particle */
+ if (gparts[k].type == swift_type_dark_matter_background &&
+ gparts[k].time_bin != time_bin_not_created) {
+
+ /* Check that it's not linked */
+ if (gparts[k].id_or_neg_offset <= 0)
+ error("Background DM gpart particle linked to something !");
+ }
+
/* We have a gas particle */
else if (gparts[k].type == swift_type_gas) {
diff --git a/src/part.h b/src/part.h
index 820834d652cf53b3dce06192d7c552ec7681f2c7..440e5e1eb3415e65d7d9cc8b611bb98de107a1b2 100644
--- a/src/part.h
+++ b/src/part.h
@@ -95,6 +95,8 @@
#include "./gravity/Default/gravity_part.h"
#elif defined(POTENTIAL_GRAVITY)
#include "./gravity/Potential/gravity_part.h"
+#elif defined(MULTI_SOFTENING_GRAVITY)
+#include "./gravity/MultiSoftening/gravity_part.h"
#else
#error "Invalid choice of gravity variant"
#endif
diff --git a/src/part_type.c b/src/part_type.c
index 1f96d4ef1db4b35a92d133e91498ea10ce472c70..5a2eee2ee2e4296f4abdd1dd3ae4927748484b15 100644
--- a/src/part_type.c
+++ b/src/part_type.c
@@ -20,5 +20,5 @@
/* This object's header. */
#include "part_type.h"
-const char* part_type_names[swift_type_count] = {"Gas", "DM", "Dummy",
- "Dummy", "Stars", "BH"};
+const char* part_type_names[swift_type_count] = {
+ "Gas", "DM", "DM_background", "Dummy", "Stars", "BH"};
diff --git a/src/part_type.h b/src/part_type.h
index 901f47193fa0e72b362c8dce5199a1d0a20526c9..8311f5d17a7b838b9577a029387a5f01fb1a0801 100644
--- a/src/part_type.h
+++ b/src/part_type.h
@@ -27,6 +27,7 @@
enum part_type {
swift_type_gas = 0,
swift_type_dark_matter = 1,
+ swift_type_dark_matter_background = 2,
swift_type_stars = 4,
swift_type_black_hole = 5,
swift_type_count
diff --git a/src/runner.c b/src/runner.c
index 22fdfaee1d422a935be5d6a5329f67d8b9c41efb..db7e512873b51a7329e19e75a763b69521efb0eb 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -2776,7 +2776,9 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gp = &gparts[k];
/* If the g-particle has no counterpart and needs to be kicked */
- if (gp->type == swift_type_dark_matter && gpart_is_starting(gp, e)) {
+ if ((gp->type == swift_type_dark_matter ||
+ gp->type == swift_type_dark_matter_background) &&
+ gpart_is_starting(gp, e)) {
const integertime_t ti_step = get_integer_timestep(gp->time_bin);
const integertime_t ti_begin =
@@ -2970,7 +2972,9 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gp = &gparts[k];
/* If the g-particle has no counterpart and needs to be kicked */
- if (gp->type == swift_type_dark_matter && gpart_is_active(gp, e)) {
+ if ((gp->type == swift_type_dark_matter ||
+ gp->type == swift_type_dark_matter_background) &&
+ gpart_is_active(gp, e)) {
const integertime_t ti_step = get_integer_timestep(gp->time_bin);
const integertime_t ti_begin =
@@ -3190,7 +3194,8 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gp = &gparts[k];
/* If the g-particle has no counterpart */
- if (gp->type == swift_type_dark_matter) {
+ if (gp->type == swift_type_dark_matter ||
+ gp->type == swift_type_dark_matter_background) {
/* need to be updated ? */
if (gpart_is_active(gp, e)) {
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index c6885746a29fd7b6bd828496316f8dad01c1b7da..b4ee8225a7aada8cf595ae7bca251d61b5226f64 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -182,12 +182,7 @@ static INLINE void runner_dopair_grav_pp_full(
const float x_i = ci_cache->x[pid];
const float y_i = ci_cache->y[pid];
const float z_i = ci_cache->z[pid];
-
- /* Some powers of the softening length */
const float h_i = ci_cache->epsilon[pid];
- const float h2_i = h_i * h_i;
- const float h_inv_i = 1.f / h_i;
- const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
/* Local accumulators for the acceleration and potential */
float a_x = 0.f, a_y = 0.f, a_z = 0.f, pot = 0.f;
@@ -197,6 +192,7 @@ static INLINE void runner_dopair_grav_pp_full(
swift_align_information(float, cj_cache->y, SWIFT_CACHE_ALIGNMENT);
swift_align_information(float, cj_cache->z, SWIFT_CACHE_ALIGNMENT);
swift_align_information(float, cj_cache->m, SWIFT_CACHE_ALIGNMENT);
+ swift_align_information(float, cj_cache->epsilon, SWIFT_CACHE_ALIGNMENT);
swift_assume_size(gcount_padded_j, VEC_SIZE);
/* Loop over every particle in the other cell. */
@@ -207,6 +203,7 @@ static INLINE void runner_dopair_grav_pp_full(
const float y_j = cj_cache->y[pjd];
const float z_j = cj_cache->z[pjd];
const float mass_j = cj_cache->m[pjd];
+ const float h_j = cj_cache->epsilon[pjd];
/* Compute the pairwise distance. */
float dx = x_j - x_i;
@@ -222,8 +219,14 @@ static INLINE void runner_dopair_grav_pp_full(
const float r2 = dx * dx + dy * dy + dz * dz;
+ /* Pick the maximal softening length of i and j */
+ const float h = max(h_i, h_j);
+ const float h2 = h * h;
+ const float h_inv = 1.f / h;
+ const float h_inv_3 = h_inv * h_inv * h_inv;
+
#ifdef SWIFT_DEBUG_CHECKS
- if (r2 == 0.f && h2_i == 0.)
+ if (r2 == 0.f && h2 == 0.)
error("Interacting particles with 0 distance and 0 softening.");
/* Check that particles have been drifted to the current time */
@@ -249,8 +252,7 @@ static INLINE void runner_dopair_grav_pp_full(
/* Interact! */
float f_ij, pot_ij;
- runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, mass_j, &f_ij,
- &pot_ij);
+ runner_iact_grav_pp_full(r2, h2, h_inv, h_inv_3, mass_j, &f_ij, &pot_ij);
/* Store it back */
a_x += f_ij * dx;
@@ -325,12 +327,7 @@ static INLINE void runner_dopair_grav_pp_truncated(
const float x_i = ci_cache->x[pid];
const float y_i = ci_cache->y[pid];
const float z_i = ci_cache->z[pid];
-
- /* Some powers of the softening length */
const float h_i = ci_cache->epsilon[pid];
- const float h2_i = h_i * h_i;
- const float h_inv_i = 1.f / h_i;
- const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
/* Local accumulators for the acceleration and potential */
float a_x = 0.f, a_y = 0.f, a_z = 0.f, pot = 0.f;
@@ -340,6 +337,7 @@ static INLINE void runner_dopair_grav_pp_truncated(
swift_align_information(float, cj_cache->y, SWIFT_CACHE_ALIGNMENT);
swift_align_information(float, cj_cache->z, SWIFT_CACHE_ALIGNMENT);
swift_align_information(float, cj_cache->m, SWIFT_CACHE_ALIGNMENT);
+ swift_align_information(float, cj_cache->epsilon, SWIFT_CACHE_ALIGNMENT);
swift_assume_size(gcount_padded_j, VEC_SIZE);
/* Loop over every particle in the other cell. */
@@ -350,6 +348,7 @@ static INLINE void runner_dopair_grav_pp_truncated(
const float y_j = cj_cache->y[pjd];
const float z_j = cj_cache->z[pjd];
const float mass_j = cj_cache->m[pjd];
+ const float h_j = cj_cache->epsilon[pjd];
/* Compute the pairwise distance. */
float dx = x_j - x_i;
@@ -363,8 +362,14 @@ static INLINE void runner_dopair_grav_pp_truncated(
const float r2 = dx * dx + dy * dy + dz * dz;
+ /* Pick the maximal softening length of i and j */
+ const float h = max(h_i, h_j);
+ const float h2 = h * h;
+ const float h_inv = 1.f / h;
+ const float h_inv_3 = h_inv * h_inv * h_inv;
+
#ifdef SWIFT_DEBUG_CHECKS
- if (r2 == 0.f && h2_i == 0.)
+ if (r2 == 0.f && h2 == 0.)
error("Interacting particles with 0 distance and 0 softening.");
/* Check that particles have been drifted to the current time */
@@ -390,8 +395,8 @@ static INLINE void runner_dopair_grav_pp_truncated(
/* Interact! */
float f_ij, pot_ij;
- runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
- r_s_inv, &f_ij, &pot_ij);
+ runner_iact_grav_pp_truncated(r2, h2, h_inv, h_inv_3, mass_j, r_s_inv,
+ &f_ij, &pot_ij);
/* Store it back */
a_x += f_ij * dx;
@@ -513,12 +518,12 @@ static INLINE void runner_dopair_grav_pm_full(
const float r_max2 = r_max_j * r_max_j;
const float theta_crit2 = e->gravity_properties->theta_crit2;
- /* Note: 1.1 to avoid FP rounding false-positives */
- if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.1, r2))
+ /* Note: 0.99 and 1.1 to avoid FP rounding false-positives */
+ if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.1, r2, 0.99 * h_i))
error(
"use_mpole[i] set when M2P accept fails CoM=[%e %e %e] pos=[%e %e "
- "%e], rmax=%e",
- CoM_j[0], CoM_j[1], CoM_j[2], x_i, y_i, z_i, r_max_j);
+ "%e], rmax=%e r=%e epsilon=%e",
+ CoM_j[0], CoM_j[1], CoM_j[2], x_i, y_i, z_i, r_max_j, sqrtf(r2), h_i);
#endif
/* Interact! */
@@ -644,8 +649,8 @@ static INLINE void runner_dopair_grav_pm_truncated(
const float r_max2 = r_max_j * r_max_j;
const float theta_crit2 = e->gravity_properties->theta_crit2;
- /* 1.1 to avoid FP rounding false-positives */
- if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.1, r2))
+ /* 0.99 and 1.1 to avoid FP rounding false-positives */
+ if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.1, r2, 0.99 * h_i))
error(
"use_mpole[i] set when M2P accept fails CoM=[%e %e %e] pos=[%e %e "
"%e], rmax=%e",
@@ -916,12 +921,7 @@ static INLINE void runner_doself_grav_pp_full(
const float x_i = ci_cache->x[pid];
const float y_i = ci_cache->y[pid];
const float z_i = ci_cache->z[pid];
-
- /* Some powers of the softening length */
const float h_i = ci_cache->epsilon[pid];
- const float h2_i = h_i * h_i;
- const float h_inv_i = 1.f / h_i;
- const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
/* Local accumulators for the acceleration */
float a_x = 0.f, a_y = 0.f, a_z = 0.f, pot = 0.f;
@@ -931,6 +931,7 @@ static INLINE void runner_doself_grav_pp_full(
swift_align_information(float, ci_cache->y, SWIFT_CACHE_ALIGNMENT);
swift_align_information(float, ci_cache->z, SWIFT_CACHE_ALIGNMENT);
swift_align_information(float, ci_cache->m, SWIFT_CACHE_ALIGNMENT);
+ swift_align_information(float, ci_cache->epsilon, SWIFT_CACHE_ALIGNMENT);
swift_assume_size(gcount_padded, VEC_SIZE);
/* Loop over every other particle in the cell. */
@@ -944,6 +945,7 @@ static INLINE void runner_doself_grav_pp_full(
const float y_j = ci_cache->y[pjd];
const float z_j = ci_cache->z[pjd];
const float mass_j = ci_cache->m[pjd];
+ const float h_j = ci_cache->epsilon[pjd];
/* Compute the pairwise (square) distance. */
/* Note: no need for periodic wrapping inside a cell */
@@ -952,8 +954,14 @@ static INLINE void runner_doself_grav_pp_full(
const float dz = z_j - z_i;
const float r2 = dx * dx + dy * dy + dz * dz;
+ /* Pick the maximal softening length of i and j */
+ const float h = max(h_i, h_j);
+ const float h2 = h * h;
+ const float h_inv = 1.f / h;
+ const float h_inv_3 = h_inv * h_inv * h_inv;
+
#ifdef SWIFT_DEBUG_CHECKS
- if (r2 == 0.f && h2_i == 0.)
+ if (r2 == 0.f && h2 == 0.)
error("Interacting particles with 0 distance and 0 softening.");
/* Check that particles have been drifted to the current time */
@@ -978,8 +986,7 @@ static INLINE void runner_doself_grav_pp_full(
/* Interact! */
float f_ij, pot_ij;
- runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, mass_j, &f_ij,
- &pot_ij);
+ runner_iact_grav_pp_full(r2, h2, h_inv, h_inv_3, mass_j, &f_ij, &pot_ij);
/* Store it back */
a_x += f_ij * dx;
@@ -1039,12 +1046,7 @@ static INLINE void runner_doself_grav_pp_truncated(
const float x_i = ci_cache->x[pid];
const float y_i = ci_cache->y[pid];
const float z_i = ci_cache->z[pid];
-
- /* Some powers of the softening length */
const float h_i = ci_cache->epsilon[pid];
- const float h2_i = h_i * h_i;
- const float h_inv_i = 1.f / h_i;
- const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
/* Local accumulators for the acceleration and potential */
float a_x = 0.f, a_y = 0.f, a_z = 0.f, pot = 0.f;
@@ -1054,6 +1056,7 @@ static INLINE void runner_doself_grav_pp_truncated(
swift_align_information(float, ci_cache->y, SWIFT_CACHE_ALIGNMENT);
swift_align_information(float, ci_cache->z, SWIFT_CACHE_ALIGNMENT);
swift_align_information(float, ci_cache->m, SWIFT_CACHE_ALIGNMENT);
+ swift_align_information(float, ci_cache->epsilon, SWIFT_CACHE_ALIGNMENT);
swift_assume_size(gcount_padded, VEC_SIZE);
/* Loop over every other particle in the cell. */
@@ -1067,6 +1070,7 @@ static INLINE void runner_doself_grav_pp_truncated(
const float y_j = ci_cache->y[pjd];
const float z_j = ci_cache->z[pjd];
const float mass_j = ci_cache->m[pjd];
+ const float h_j = ci_cache->epsilon[pjd];
/* Compute the pairwise (square) distance. */
/* Note: no need for periodic wrapping inside a cell */
@@ -1076,8 +1080,14 @@ static INLINE void runner_doself_grav_pp_truncated(
const float r2 = dx * dx + dy * dy + dz * dz;
+ /* Pick the maximal softening length of i and j */
+ const float h = max(h_i, h_j);
+ const float h2 = h * h;
+ const float h_inv = 1.f / h;
+ const float h_inv_3 = h_inv * h_inv * h_inv;
+
#ifdef SWIFT_DEBUG_CHECKS
- if (r2 == 0.f && h2_i == 0.)
+ if (r2 == 0.f && h2 == 0.)
error("Interacting particles with 0 distance and 0 softening.");
/* Check that particles have been drifted to the current time */
@@ -1102,8 +1112,8 @@ static INLINE void runner_doself_grav_pp_truncated(
/* Interact! */
float f_ij, pot_ij;
- runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
- r_s_inv, &f_ij, &pot_ij);
+ runner_iact_grav_pp_truncated(r2, h2, h_inv, h_inv_3, mass_j, r_s_inv,
+ &f_ij, &pot_ij);
/* Store it back */
a_x += f_ij * dx;
@@ -1575,7 +1585,9 @@ static INLINE void runner_dopair_recursive_grav(struct runner *r,
* option... */
/* Can we use M-M interactions ? */
- if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2)) {
+ if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2,
+ multi_i->m_pole.max_softening,
+ multi_j->m_pole.max_softening)) {
/* Go M-M */
runner_dopair_grav_mm(r, ci, cj);
@@ -1797,7 +1809,9 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
/* Are we in charge of this cell pair? */
if (gravity_M2L_accept(multi_top->r_max_rebuild, multi_j->r_max_rebuild,
- theta_crit2, r2_rebuild)) {
+ theta_crit2, r2_rebuild,
+ multi_top->m_pole.max_softening,
+ multi_j->m_pole.max_softening)) {
/* Call the PM interaction fucntion on the active sub-cells of ci */
runner_dopair_grav_mm_nonsym(r, ci, cj);
diff --git a/src/serial_io.c b/src/serial_io.c
index 0d367b1cf423d0471a0f6235fa3a1c6877e9f279..ca832a4beb16f8f39760c3c0a2950ecd003cc02b 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -460,6 +460,8 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
* @param bparts (output) Array of #bpart particles.
* @param Ngas (output) The number of #part read from the file on that node.
* @param Ngparts (output) The number of #gpart read from the file on that node.
+ * @param Ngparts_background (output) The number of background #gpart (type 2)
+ * read from the file on that node.
* @param Nstars (output) The number of #spart read from the file on that node.
* @param Nblackholes (output) The number of #bpart read from the file on that
* node.
@@ -492,12 +494,12 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
void read_ic_serial(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, struct bpart** bparts, size_t* Ngas,
- size_t* Ngparts, size_t* Nstars, size_t* Nblackholes,
- int* flag_entropy, int with_hydro, int with_gravity,
- int with_stars, int with_black_holes, int cleanup_h,
- int cleanup_sqrt_a, double h, double a, int mpi_rank,
- int mpi_size, MPI_Comm comm, MPI_Info info, int n_threads,
- int dry_run) {
+ size_t* Ngparts, size_t* Ngparts_background, size_t* Nstars,
+ size_t* Nblackholes, int* flag_entropy, int with_hydro,
+ int with_gravity, int with_stars, int with_black_holes,
+ int cleanup_h, int cleanup_sqrt_a, double h, double a,
+ int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
+ int n_threads, int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
@@ -510,11 +512,13 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
long long offset[swift_type_count] = {0};
int dimension = 3; /* Assume 3D if nothing is specified */
size_t Ndm = 0;
+ size_t Ndm_background = 0;
struct unit_system* ic_units =
(struct unit_system*)malloc(sizeof(struct unit_system));
/* Initialise counters */
- *Ngas = 0, *Ngparts = 0, *Nstars = 0, *Nblackholes = 0;
+ *Ngas = 0, *Ngparts = 0, *Ngparts_background = 0, *Nstars = 0,
+ *Nblackholes = 0;
/* First read some information about the content */
if (mpi_rank == 0) {
@@ -675,11 +679,14 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
/* Allocate memory to store all gravity particles */
if (with_gravity) {
- Ndm = N[1];
+ Ndm = N[swift_type_dark_matter];
+ Ndm_background = N[swift_type_dark_matter_background];
*Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
N[swift_type_dark_matter] +
+ N[swift_type_dark_matter_background] +
(with_stars ? N[swift_type_stars] : 0) +
(with_black_holes ? N[swift_type_black_hole] : 0);
+ *Ngparts_background = Ndm_background;
if (swift_memalign("gparts", (void**)gparts, gpart_align,
*Ngparts * sizeof(struct gpart)) != 0)
error("Error while allocating memory for gravity particles");
@@ -740,6 +747,13 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
}
break;
+ case swift_type_dark_matter_background:
+ if (with_gravity) {
+ Nparticles = Ndm_background;
+ darkmatter_read_particles(*gparts + Ndm, list, &num_fields);
+ }
+ break;
+
case swift_type_stars:
if (with_stars) {
Nparticles = *Nstars;
@@ -789,17 +803,23 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
/* Prepare the DM particles */
io_prepare_dm_gparts(&tp, *gparts, Ndm);
+ /* Prepare the DM background particles */
+ io_prepare_dm_background_gparts(&tp, *gparts + Ndm, Ndm_background);
+
/* Duplicate the hydro particles into gparts */
- if (with_hydro) io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas, Ndm);
+ if (with_hydro)
+ io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas,
+ Ndm + Ndm_background);
/* Duplicate the stars particles into gparts */
if (with_stars)
- io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
+ io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars,
+ Ndm + Ndm_background + *Ngas);
/* Duplicate the black holes particles into gparts */
if (with_black_holes)
io_duplicate_black_holes_gparts(&tp, *bparts, *gparts, *Nblackholes,
- Ndm + *Ngas + *Nstars);
+ Ndm + Ndm_background + *Ngas + *Nstars);
threadpool_clean(&tp);
}
@@ -847,6 +867,7 @@ void write_output_serial(struct engine* e, const char* baseName,
const int with_cooling = e->policy & engine_policy_cooling;
const int with_temperature = e->policy & engine_policy_temperature;
const int with_fof = e->policy & engine_policy_fof;
+ const int with_DM_background = e->s->with_DM_background;
#ifdef HAVE_VELOCIRAPTOR
const int with_stf = (e->policy & engine_policy_structure_finding) &&
(e->s->gpart_group_data != NULL);
@@ -864,7 +885,13 @@ void write_output_serial(struct engine* e, const char* baseName,
// const size_t Nbaryons = Ngas + Nstars;
// const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
- /* Number of particles that we will write */
+ size_t Ndm_background = 0;
+ if (with_DM_background) {
+ Ndm_background = io_count_dm_background_gparts(gparts, Ntot);
+ }
+
+ /* Number of particles that we will write
+ * Recall that background particles are never inhibited and have no extras */
const size_t Ntot_written =
e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
const size_t Ngas_written =
@@ -876,7 +903,7 @@ void write_output_serial(struct engine* e, const char* baseName,
const size_t Nbaryons_written =
Ngas_written + Nstars_written + Nblackholes_written;
const size_t Ndm_written =
- Ntot_written > 0 ? Ntot_written - Nbaryons_written : 0;
+ Ntot_written > 0 ? Ntot_written - Nbaryons_written - Ndm_background : 0;
/* File name */
char fileName[FILENAME_BUFFER_SIZE];
@@ -888,7 +915,8 @@ void write_output_serial(struct engine* e, const char* baseName,
e->snapshot_output_count);
/* Compute offset in the file and total number of particles */
- size_t N[swift_type_count] = {Ngas_written, Ndm_written, 0, 0,
+ size_t N[swift_type_count] = {Ngas_written, Ndm_written,
+ Ndm_background, 0,
Nstars_written, Nblackholes_written};
long long N_total[swift_type_count] = {0};
long long offset[swift_type_count] = {0};
@@ -1242,7 +1270,8 @@ void write_output_serial(struct engine* e, const char* baseName,
case swift_type_dark_matter: {
if (Ntot == Ndm_written) {
- /* This is a DM-only run without inhibited particles */
+ /* This is a DM-only run without background or inhibited particles
+ */
Nparticles = Ntot;
darkmatter_write_particles(gparts, list, &num_fields);
if (with_fof) {
@@ -1293,6 +1322,45 @@ void write_output_serial(struct engine* e, const char* baseName,
}
} break;
+ case swift_type_dark_matter_background: {
+
+ /* Ok, we need to fish out the particles we want */
+ Nparticles = Ndm_background;
+
+ /* Allocate temporary array */
+ if (swift_memalign("gparts_written", (void**)&gparts_written,
+ gpart_align,
+ Ndm_background * sizeof(struct gpart)) != 0)
+ error("Error while allocating temporart memory for gparts");
+
+ if (with_stf) {
+ if (swift_memalign("gpart_group_written",
+ (void**)&gpart_group_data_written, gpart_align,
+ Ndm_background *
+ sizeof(struct velociraptor_gpart_data)) !=
+ 0)
+ error(
+ "Error while allocating temporart memory for gparts STF "
+ "data");
+ }
+
+ /* Collect the non-inhibited DM particles from gpart */
+ io_collect_gparts_background_to_write(
+ gparts, e->s->gpart_group_data, gparts_written,
+ gpart_group_data_written, Ntot, Ndm_background, with_stf);
+
+ /* Select the fields to write */
+ darkmatter_write_particles(gparts_written, list, &num_fields);
+ if (with_fof) {
+ num_fields += fof_write_gparts(gparts_written, list + num_fields);
+ }
+ if (with_stf) {
+ num_fields += velociraptor_write_gparts(gpart_group_data_written,
+ list + num_fields);
+ }
+
+ } break;
+
case swift_type_stars: {
if (Nstars == Nstars_written) {
diff --git a/src/serial_io.h b/src/serial_io.h
index 7d994181b135e9758a9e3a52971366c2d871f627..c6ca7ad66df038791fbd62277f406a7f34340e04 100644
--- a/src/serial_io.h
+++ b/src/serial_io.h
@@ -38,12 +38,12 @@
void read_ic_serial(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, struct bpart** bparts, size_t* Ngas,
- size_t* Ngparts, size_t* Nstars, size_t* Nblackholes,
- int* flag_entropy, int with_hydro, int with_gravity,
- int with_stars, int with_black_holes, int cleanup_h,
- int cleanup_sqrt_a, double h, double a, int mpi_rank,
- int mpi_size, MPI_Comm comm, MPI_Info info, int n_threads,
- int dry_run);
+ size_t* Ngparts, size_t* Ngparts_background, size_t* Nstars,
+ size_t* Nblackholes, int* flag_entropy, int with_hydro,
+ int with_gravity, int with_stars, int with_black_holes,
+ int cleanup_h, int cleanup_sqrt_a, double h, double a,
+ int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
+ int n_threads, int dry_run);
void write_output_serial(struct engine* e, const char* baseName,
const struct unit_system* internal_units,
@@ -56,6 +56,7 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
int mpi_rank, long long offset,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units);
-#endif
+
+#endif /* HAVE_HDF5 && WITH_MPI && !HAVE_PARALLEL_HDF5 */
#endif /* SWIFT_SERIAL_IO_H */
diff --git a/src/single_io.c b/src/single_io.c
index 76a5630120df5e946093e9f046503f73bce706b6..017f9498ffc7d11e1a03d27eca733d2d417b9777 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -404,11 +404,11 @@ void read_ic_single(const char* fileName,
const struct unit_system* internal_units, double dim[3],
struct part** parts, struct gpart** gparts,
struct spart** sparts, struct bpart** bparts, size_t* Ngas,
- size_t* Ngparts, size_t* Nstars, size_t* Nblackholes,
- int* flag_entropy, int with_hydro, int with_gravity,
- int with_stars, int with_black_holes, int cleanup_h,
- int cleanup_sqrt_a, double h, double a, int n_threads,
- int dry_run) {
+ size_t* Ngparts, size_t* Ngparts_background, size_t* Nstars,
+ size_t* Nblackholes, int* flag_entropy, int with_hydro,
+ int with_gravity, int with_stars, int with_black_holes,
+ int cleanup_h, int cleanup_sqrt_a, double h, double a,
+ int n_threads, int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
@@ -419,9 +419,11 @@ void read_ic_single(const char* fileName,
size_t N[swift_type_count] = {0};
int dimension = 3; /* Assume 3D if nothing is specified */
size_t Ndm = 0;
+ size_t Ndm_background = 0;
/* Initialise counters */
- *Ngas = 0, *Ngparts = 0, *Nstars = 0, *Nblackholes = 0;
+ *Ngas = 0, *Ngparts = 0, *Ngparts_background = 0, *Nstars = 0,
+ *Nblackholes = 0;
/* Open file */
/* message("Opening file '%s' as IC.", fileName); */
@@ -563,10 +565,13 @@ void read_ic_single(const char* fileName,
/* Allocate memory to store all gravity particles */
if (with_gravity) {
Ndm = N[swift_type_dark_matter];
+ Ndm_background = N[swift_type_dark_matter_background];
*Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
N[swift_type_dark_matter] +
+ N[swift_type_dark_matter_background] +
(with_stars ? N[swift_type_stars] : 0) +
(with_black_holes ? N[swift_type_black_hole] : 0);
+ *Ngparts_background = Ndm_background;
if (swift_memalign("gparts", (void**)gparts, gpart_align,
*Ngparts * sizeof(struct gpart)) != 0)
error("Error while allocating memory for gravity particles");
@@ -615,6 +620,13 @@ void read_ic_single(const char* fileName,
}
break;
+ case swift_type_dark_matter_background:
+ if (with_gravity) {
+ Nparticles = Ndm_background;
+ darkmatter_read_particles(*gparts + Ndm, list, &num_fields);
+ }
+ break;
+
case swift_type_stars:
if (with_stars) {
Nparticles = *Nstars;
@@ -653,17 +665,23 @@ void read_ic_single(const char* fileName,
/* Prepare the DM particles */
io_prepare_dm_gparts(&tp, *gparts, Ndm);
+ /* Prepare the DM background particles */
+ io_prepare_dm_background_gparts(&tp, *gparts + Ndm, Ndm_background);
+
/* Duplicate the hydro particles into gparts */
- if (with_hydro) io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas, Ndm);
+ if (with_hydro)
+ io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas,
+ Ndm + Ndm_background);
/* Duplicate the star particles into gparts */
if (with_stars)
- io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
+ io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars,
+ Ndm + Ndm_background + *Ngas);
/* Duplicate the black hole particles into gparts */
if (with_black_holes)
io_duplicate_black_holes_gparts(&tp, *bparts, *gparts, *Nblackholes,
- Ndm + *Ngas + *Nstars);
+ Ndm + Ndm_background + *Ngas + *Nstars);
threadpool_clean(&tp);
}
@@ -709,6 +727,7 @@ void write_output_single(struct engine* e, const char* baseName,
const int with_cooling = e->policy & engine_policy_cooling;
const int with_temperature = e->policy & engine_policy_temperature;
const int with_fof = e->policy & engine_policy_fof;
+ const int with_DM_background = e->s->with_DM_background;
#ifdef HAVE_VELOCIRAPTOR
const int with_stf = (e->policy & engine_policy_structure_finding) &&
(e->s->gpart_group_data != NULL);
@@ -724,7 +743,13 @@ void write_output_single(struct engine* e, const char* baseName,
// const size_t Nbaryons = Ngas + Nstars;
// const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
- /* Number of particles that we will write */
+ size_t Ndm_background = 0;
+ if (with_DM_background) {
+ Ndm_background = io_count_dm_background_gparts(gparts, Ntot);
+ }
+
+ /* Number of particles that we will write
+ * Recall that background particles are never inhibited and have no extras */
const size_t Ntot_written =
e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
const size_t Ngas_written =
@@ -736,11 +761,12 @@ void write_output_single(struct engine* e, const char* baseName,
const size_t Nbaryons_written =
Ngas_written + Nstars_written + Nblackholes_written;
const size_t Ndm_written =
- Ntot_written > 0 ? Ntot_written - Nbaryons_written : 0;
+ Ntot_written > 0 ? Ntot_written - Nbaryons_written - Ndm_background : 0;
/* Format things in a Gadget-friendly array */
long long N_total[swift_type_count] = {
- (long long)Ngas_written, (long long)Ndm_written, 0, 0,
+ (long long)Ngas_written, (long long)Ndm_written,
+ (long long)Ndm_background, 0,
(long long)Nstars_written, (long long)Nblackholes_written};
/* File name */
@@ -1040,7 +1066,7 @@ void write_output_single(struct engine* e, const char* baseName,
case swift_type_dark_matter: {
if (Ntot == Ndm_written) {
- /* This is a DM-only run without inhibited particles */
+ /* This is a DM-only run without background or inhibited particles */
N = Ntot;
darkmatter_write_particles(gparts, list, &num_fields);
if (with_fof) {
@@ -1088,6 +1114,43 @@ void write_output_single(struct engine* e, const char* baseName,
}
} break;
+ case swift_type_dark_matter_background: {
+
+ /* Ok, we need to fish out the particles we want */
+ N = Ndm_background;
+
+ /* Allocate temporary array */
+ if (swift_memalign("gparts_written", (void**)&gparts_written,
+ gpart_align,
+ Ndm_background * sizeof(struct gpart)) != 0)
+ error("Error while allocating temporart memory for gparts");
+
+ if (with_stf) {
+ if (swift_memalign(
+ "gpart_group_written", (void**)&gpart_group_data_written,
+ gpart_align,
+ Ndm_background * sizeof(struct velociraptor_gpart_data)) != 0)
+ error(
+ "Error while allocating temporart memory for gparts STF "
+ "data");
+ }
+
+ /* Collect the non-inhibited DM particles from gpart */
+ io_collect_gparts_background_to_write(
+ gparts, e->s->gpart_group_data, gparts_written,
+ gpart_group_data_written, Ntot, Ndm_background, with_stf);
+
+ /* Select the fields to write */
+ darkmatter_write_particles(gparts_written, list, &num_fields);
+ if (with_fof) {
+ num_fields += fof_write_gparts(gparts_written, list + num_fields);
+ }
+ if (with_stf) {
+ num_fields += velociraptor_write_gparts(gpart_group_data_written,
+ list + num_fields);
+ }
+ } break;
+
case swift_type_stars: {
if (Nstars == Nstars_written) {
@@ -1224,4 +1287,4 @@ void write_output_single(struct engine* e, const char* baseName,
e->snapshot_output_count++;
}
-#endif /* HAVE_HDF5 */
+#endif /* HAVE_HDF5 && !WITH_MPI */
diff --git a/src/single_io.h b/src/single_io.h
index 9ff04c893378d7f8c01621e7dbf387dc939d0157..47b575272c2fbe5b170eb4370743cd6830b99ae0 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -34,11 +34,11 @@ void read_ic_single(const char* fileName,
const struct unit_system* internal_units, double dim[3],
struct part** parts, struct gpart** gparts,
struct spart** sparts, struct bpart** bparts, size_t* Ngas,
- size_t* Ndm, size_t* Nstars, size_t* Nblackholes,
- int* flag_entropy, int with_hydro, int with_gravity,
- int with_stars, int with_black_holes, int cleanup_h,
- int cleanup_sqrt_a, double h, double a, int nr_threads,
- int dry_run);
+ size_t* Ndm, size_t* Ndm_background, size_t* Nstars,
+ size_t* Nblackholes, int* flag_entropy, int with_hydro,
+ int with_gravity, int with_stars, int with_black_holes,
+ int cleanup_h, int cleanup_sqrt_a, double h, double a,
+ int nr_threads, int dry_run);
void write_output_single(struct engine* e, const char* baseName,
const struct unit_system* internal_units,
diff --git a/src/space.c b/src/space.c
index 6594fab935d55e969faa7482e71f4a2b73954399..522347d03abff370b7e175f322166395d38e1321 100644
--- a/src/space.c
+++ b/src/space.c
@@ -1931,7 +1931,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* Check that the multipole construction went OK */
if (s->with_self_gravity)
for (int k = 0; k < s->nr_cells; k++)
- cell_check_multipole(&s->cells_top[k]);
+ cell_check_multipole(&s->cells_top[k], s->e->gravity_properties);
#endif
/* Clean up any stray sort indices in the cell buffer. */
@@ -3606,7 +3606,8 @@ void space_split_recursive(struct space *s, struct cell *c,
if (s->with_self_gravity) {
if (gcount > 0) {
- gravity_P2M(c->grav.multipole, c->grav.parts, c->grav.count);
+ gravity_P2M(c->grav.multipole, c->grav.parts, c->grav.count,
+ e->gravity_properties);
} else {
@@ -3948,6 +3949,9 @@ void space_synchronize_particle_positions_mapper(void *map_data, int nr_gparts,
if (gp->type == swift_type_dark_matter)
continue;
+ else if (gp->type == swift_type_dark_matter_background)
+ continue;
+
else if (gp->type == swift_type_gas) {
/* Get its gassy friend */
@@ -3991,6 +3995,10 @@ void space_synchronize_particle_positions_mapper(void *map_data, int nr_gparts,
gp->mass = bp->mass;
}
+
+ else {
+ error("Invalid type!");
+ }
}
}
@@ -4526,6 +4534,7 @@ void space_convert_quantities(struct space *s, int verbose) {
* @param hydro flag whether we are doing hydro or not?
* @param self_gravity flag whether we are doing gravity or not?
* @param star_formation flag whether we are doing star formation or not?
+ * @param DM_background Are we running with some DM background particles?
* @param verbose Print messages to stdout or not.
* @param dry_run If 1, just initialise stuff, don't do anything with the parts.
*
@@ -4540,7 +4549,8 @@ void space_init(struct space *s, struct swift_params *params,
struct bpart *bparts, size_t Npart, size_t Ngpart,
size_t Nspart, size_t Nbpart, int periodic, int replicate,
int generate_gas_in_ics, int hydro, int self_gravity,
- int star_formation, int verbose, int dry_run) {
+ int star_formation, int DM_background, int verbose,
+ int dry_run) {
/* Clean-up everything */
bzero(s, sizeof(struct space));
@@ -4553,6 +4563,7 @@ void space_init(struct space *s, struct swift_params *params,
s->with_self_gravity = self_gravity;
s->with_hydro = hydro;
s->with_star_formation = star_formation;
+ s->with_DM_background = DM_background;
s->nr_parts = Npart;
s->nr_gparts = Ngpart;
s->nr_sparts = Nspart;
@@ -4603,6 +4614,9 @@ void space_init(struct space *s, struct swift_params *params,
error("Value of 'InitialConditions:replicate' (%d) is too small",
replicate);
if (replicate > 1) {
+ if (DM_background)
+ error("Can't replicate the space if background DM particles are in use.");
+
space_replicate(s, replicate, verbose);
parts = s->parts;
gparts = s->gparts;
@@ -4841,16 +4855,19 @@ void space_replicate(struct space *s, int replicate, int verbose) {
const size_t nr_parts = s->nr_parts;
const size_t nr_gparts = s->nr_gparts;
const size_t nr_sparts = s->nr_sparts;
- const size_t nr_dm = nr_gparts - nr_parts - nr_sparts;
+ const size_t nr_bparts = s->nr_bparts;
+ const size_t nr_dm = nr_gparts - nr_parts - nr_sparts - nr_bparts;
s->size_parts = s->nr_parts = nr_parts * factor;
s->size_gparts = s->nr_gparts = nr_gparts * factor;
s->size_sparts = s->nr_sparts = nr_sparts * factor;
+ s->size_bparts = s->nr_bparts = nr_bparts * factor;
/* Allocate space for new particles */
struct part *parts = NULL;
struct gpart *gparts = NULL;
struct spart *sparts = NULL;
+ struct bpart *bparts = NULL;
if (swift_memalign("parts", (void **)&parts, part_align,
s->nr_parts * sizeof(struct part)) != 0)
@@ -4864,6 +4881,10 @@ void space_replicate(struct space *s, int replicate, int verbose) {
s->nr_sparts * sizeof(struct spart)) != 0)
error("Failed to allocate new spart array.");
+ if (swift_memalign("bparts", (void **)&bparts, bpart_align,
+ s->nr_bparts * sizeof(struct bpart)) != 0)
+ error("Failed to allocate new bpart array.");
+
/* Replicate everything */
for (int i = 0; i < replicate; ++i) {
for (int j = 0; j < replicate; ++j) {
@@ -4875,6 +4896,8 @@ void space_replicate(struct space *s, int replicate, int verbose) {
nr_parts * sizeof(struct part));
memcpy(sparts + offset * nr_sparts, s->sparts,
nr_sparts * sizeof(struct spart));
+ memcpy(bparts + offset * nr_bparts, s->bparts,
+ nr_bparts * sizeof(struct bpart));
memcpy(gparts + offset * nr_gparts, s->gparts,
nr_gparts * sizeof(struct gpart));
@@ -4896,6 +4919,11 @@ void space_replicate(struct space *s, int replicate, int verbose) {
sparts[n].x[1] += shift[1];
sparts[n].x[2] += shift[2];
}
+ for (size_t n = offset * nr_bparts; n < (offset + 1) * nr_bparts; ++n) {
+ bparts[n].x[0] += shift[0];
+ bparts[n].x[1] += shift[1];
+ bparts[n].x[2] += shift[2];
+ }
/* Set the correct links (recall gpart are sorted by type at start-up):
first DM (unassociated gpart), then gas, then stars */
@@ -4917,6 +4945,16 @@ void space_replicate(struct space *s, int replicate, int verbose) {
gparts[offset_gpart + n].id_or_neg_offset = -(offset_spart + n);
}
}
+ if (nr_bparts > 0 && nr_gparts > 0) {
+ const size_t offset_bpart = offset * nr_bparts;
+ const size_t offset_gpart =
+ offset * nr_gparts + nr_dm + nr_parts + nr_sparts;
+
+ for (size_t n = 0; n < nr_bparts; ++n) {
+ bparts[offset_bpart + n].gpart = &gparts[offset_gpart + n];
+ gparts[offset_gpart + n].id_or_neg_offset = -(offset_bpart + n);
+ }
+ }
}
}
}
@@ -4925,9 +4963,11 @@ void space_replicate(struct space *s, int replicate, int verbose) {
swift_free("parts", s->parts);
swift_free("gparts", s->gparts);
swift_free("sparts", s->sparts);
+ swift_free("bparts", s->bparts);
s->parts = parts;
s->gparts = gparts;
s->sparts = sparts;
+ s->bparts = bparts;
/* Finally, update the domain size */
s->dim[0] *= replicate;
diff --git a/src/space.h b/src/space.h
index 1b31f0495890a68074225cc68c6acabd1abb53e8..0b332716645e733636b7ab0da57a0a31b28e3d31 100644
--- a/src/space.h
+++ b/src/space.h
@@ -103,6 +103,9 @@ struct space {
/*! Are we doing star formation? */
int with_star_formation;
+ /*! Are we running with some DM background particles? */
+ int with_DM_background;
+
/*! Width of the top-level cells. */
double width[3];
@@ -304,7 +307,8 @@ void space_init(struct space *s, struct swift_params *params,
struct bpart *bparts, size_t Npart, size_t Ngpart,
size_t Nspart, size_t Nbpart, int periodic, int replicate,
int generate_gas_in_ics, int hydro, int gravity,
- int star_formation, int verbose, int dry_run);
+ int star_formation, int DM_background, int verbose,
+ int dry_run);
void space_sanitize(struct space *s);
void space_map_cells_pre(struct space *s, int full,
void (*fun)(struct cell *c, void *data), void *data);
diff --git a/src/velociraptor_interface.c b/src/velociraptor_interface.c
index 47c107d047e3abfcac32571f780cf52649a9c38d..8cb65f6ac13fb6d4f813d68656363640b082123f 100644
--- a/src/velociraptor_interface.c
+++ b/src/velociraptor_interface.c
@@ -283,6 +283,13 @@ void velociraptor_convert_particles_mapper(void *map_data, int nr_gparts,
swift_parts[i].T = 0.f;
break;
+ case swift_type_dark_matter_background:
+
+ swift_parts[i].id = gparts[i].id_or_neg_offset;
+ swift_parts[i].u = 0.f;
+ swift_parts[i].T = 0.f;
+ break;
+
default:
error("Particle type not handled by VELOCIraptor.");
}
@@ -334,13 +341,19 @@ void velociraptor_init(struct engine *e) {
} else {
sim_info.icosmologicalsim = 0;
}
- sim_info.izoomsim = 0;
+
+ /* Are we running a zoom? */
+ if (e->s->with_DM_background) {
+ sim_info.izoomsim = 1;
+ } else {
+ sim_info.izoomsim = 0;
+ }
/* Tell VELOCIraptor what we have in the simulation */
sim_info.idarkmatter = (e->total_nr_gparts - e->total_nr_parts > 0);
sim_info.igas = (e->policy & engine_policy_hydro);
sim_info.istar = (e->policy & engine_policy_stars);
- sim_info.ibh = 0; // sim_info.ibh = (e->policy&engine_policy_bh);
+ sim_info.ibh = (e->policy & engine_policy_black_holes);
sim_info.iother = 0;
/* Be nice, talk! */
@@ -443,14 +456,45 @@ void velociraptor_invoke(struct engine *e, const int linked_with_snap) {
/* Are we running with cosmology? */
if (e->policy & engine_policy_cosmology) {
sim_info.icosmologicalsim = 1;
- sim_info.izoomsim = 0;
- const size_t total_nr_baryons = e->total_nr_parts + e->total_nr_sparts;
- const size_t total_nr_dmparts = e->total_nr_gparts - total_nr_baryons;
- sim_info.interparticlespacing = sim_info.period / cbrt(total_nr_dmparts);
+
+ /* Are we running a zoom? */
+ if (e->s->with_DM_background) {
+ sim_info.izoomsim = 1;
+ } else {
+ sim_info.izoomsim = 0;
+ }
+
+ /* Collect the mass of the non-background gpart */
+ double high_res_DM_mass = 0.;
+ for (size_t i = 0; i < e->s->nr_gparts; ++i) {
+ const struct gpart *gp = &e->s->gparts[i];
+ if (gp->type == swift_type_dark_matter &&
+ gp->time_bin != time_bin_inhibited &&
+ gp->time_bin != time_bin_not_created) {
+ high_res_DM_mass = gp->mass;
+ break;
+ }
+ }
+
+#ifdef WITH_MPI
+ /* We need to all-reduce this in case one of the nodes had 0 DM particles.
+ */
+ MPI_Allreduce(MPI_IN_PLACE, &high_res_DM_mass, 1, MPI_DOUBLE, MPI_MAX,
+ MPI_COMM_WORLD);
+#endif
+
+ /* Linking length based on the mean DM inter-particle separation
+ * in the zoom region and assuming the mean density of the Universe
+ * is used in the zoom region. */
+ const double mean_matter_density =
+ e->cosmology->Omega_m * e->cosmology->critical_density_0;
+ sim_info.interparticlespacing =
+ cbrt(high_res_DM_mass / mean_matter_density);
+
} else {
- sim_info.icosmologicalsim = 0;
sim_info.izoomsim = 0;
- sim_info.interparticlespacing = -1;
+ sim_info.icosmologicalsim = 0;
+ sim_info.interparticlespacing = -1.;
}
/* Set the spatial extent of the simulation volume */
diff --git a/tests/testGravityDerivatives.c b/tests/testGravityDerivatives.c
index f31967de7075bccfb2c7fb19c1ba262aa12da54f..032460d9f519c455e485e0dddfc5f8742138d60a 100644
--- a/tests/testGravityDerivatives.c
+++ b/tests/testGravityDerivatives.c
@@ -955,12 +955,11 @@ int main(int argc, char* argv[]) {
const double r_inv = 1. / sqrt(r2);
const double r = r2 * r_inv;
const double eps = r / 10.;
- const double eps_inv = 1. / eps;
/* Compute all derivatives */
struct potential_derivatives_M2L pot;
- potential_derivatives_compute_M2L(dx, dy, dz, r2, r_inv, eps, eps_inv,
- periodic, r_s_inv, &pot);
+ potential_derivatives_compute_M2L(dx, dy, dz, r2, r_inv, eps, periodic,
+ r_s_inv, &pot);
/* Minimal value we care about */
const double min = 1e-9;
diff --git a/tests/testPotentialPair.c b/tests/testPotentialPair.c
index 064c86d42f8df907d1ffaaab164b6a2f8b534b19..d5fbda36a9ef79352f79627b5cef908030401da2 100644
--- a/tests/testPotentialPair.c
+++ b/tests/testPotentialPair.c
@@ -125,7 +125,8 @@ int main(int argc, char *argv[]) {
struct gravity_props props;
props.theta_crit2 = 0.;
- props.epsilon_cur = eps;
+ props.epsilon_DM_cur = eps;
+ props.epsilon_baryon_cur = eps;
e.gravity_properties = &props;
struct runner r;
@@ -386,7 +387,7 @@ int main(int argc, char *argv[]) {
/* Now let's make a multipole out of it. */
gravity_reset(ci.grav.multipole);
- gravity_P2M(ci.grav.multipole, ci.grav.parts, ci.grav.count);
+ gravity_P2M(ci.grav.multipole, ci.grav.parts, ci.grav.count, &props);
gravity_multipole_print(&ci.grav.multipole->m_pole);
diff --git a/tests/testPotentialSelf.c b/tests/testPotentialSelf.c
index 10eb499570a591daaf0de2e011f2346077905e8e..ff55558aff77441458e822cfa11922bd2dabfa06 100644
--- a/tests/testPotentialSelf.c
+++ b/tests/testPotentialSelf.c
@@ -115,7 +115,8 @@ int main(int argc, char *argv[]) {
struct gravity_props props;
props.a_smooth = 1.25;
- props.epsilon_cur = eps;
+ props.epsilon_DM_cur = eps;
+ props.epsilon_baryon_cur = eps;
e.gravity_properties = &props;
struct runner r;
diff --git a/tests/testReading.c b/tests/testReading.c
index b2cf743a066920c3d28ea8768334a6a8b1c9b5f0..47f8d0dfdc043985c6344abecf31e605226e6899 100644
--- a/tests/testReading.c
+++ b/tests/testReading.c
@@ -28,7 +28,7 @@
int main(int argc, char *argv[]) {
- size_t Ngas = 0, Ngpart = 0, Nspart = 0, Nbpart = 0;
+ size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0;
int flag_entropy_ICs = -1;
int i, j, k;
double dim[3];
@@ -51,8 +51,15 @@ int main(int argc, char *argv[]) {
/* Read data */
read_ic_single("input.hdf5", &us, dim, &parts, &gparts, &sparts, &bparts,
- &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs, 1, 1, 0,
- 0, 0, 0, 1., 1., 1, 0);
+ &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
+ &flag_entropy_ICs,
+ /*with_hydro=*/1,
+ /*with_gravity=*/1,
+ /*with_stars=*/0,
+ /*with_black_holes=*/0,
+ /*cleanup_h=*/0,
+ /*cleanup_sqrt_a=*/0,
+ /*h=*/1., /*a=*/1., /*n_threads=*/1, /*dry_run=*/0);
/* Check global properties read are correct */
assert(dim[0] == boxSize);
diff --git a/tests/testSelectOutput.c b/tests/testSelectOutput.c
index 0be250e48cb2748093c8cd8c1381303c3060c9ba..9797406d91b74c7e9effe14da7b4f0100d901148 100644
--- a/tests/testSelectOutput.c
+++ b/tests/testSelectOutput.c
@@ -84,7 +84,7 @@ int main(int argc, char *argv[]) {
clocks_set_cpufreq(cpufreq);
char *base_name = "testSelectOutput";
- size_t Ngas = 0, Ngpart = 0, Nspart = 0, Nbpart = 0;
+ size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0;
int flag_entropy_ICs = -1;
int periodic = 1;
double dim[3];
@@ -112,8 +112,15 @@ int main(int argc, char *argv[]) {
/* Read data */
message("Reading initial conditions.");
read_ic_single("input.hdf5", &us, dim, &parts, &gparts, &sparts, &bparts,
- &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs, 1, 0, 0,
- 0, 0, 0, 1., 1., 1, 0);
+ &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
+ &flag_entropy_ICs,
+ /*with_hydro=*/1,
+ /*with_gravity=*/0,
+ /*with_stars=*/0,
+ /*with_black_holes=*/0,
+ /*cleanup_h=*/0,
+ /*cleanup_sqrt_a=*/0,
+ /*h=*/1., /*a=*/1., /*n_threads=*/1, /*dry_run=*/0);
/* pseudo initialization of the space */
message("Initialization of the space.");