Sphenix pressure floor

doc/RTD/source/HydroSchemes/sphenix_sph.rst

@@ -67,3 +67,14 @@ There is also a compile-time parameter, ``viscosity_beta`` that we set to
 ``hydro_props_default_viscosity_alpha_feedback_reset = 2.0`` and the
 diffusion is set to ``hydro_props_default_diffusion_alpha_feedback_reset =
 0.0``. These can be changed in ``src/hydro/SPHENIX/hydro_parameters.h``.
+
+Pressure Floor
+~~~~~~~~~~~~~~
+
+The pressure floor is implemented for this scheme. The 
+pressure floor is used in the sound-speed and hence the
+time-step will be lower in this case than if the pressure
+floor was not taken into account. The only other impact
+of the pressure floor is in the force loop, where the
+sub-grid floor pressure is used instead of the pressure
+calculated from the internal energy and density.

doc/RTD/source/SubgridModels/GEAR/index.rst

@@ -19,7 +19,7 @@ where :math:`\rho` is the density, :math:`\gamma` the adiabatic index, :math:`G`
 :math:`\sigma` the velocity dispersion.
 
 
-This must be directly implemented into the hydro schemes, therefore only a subset of schemes (Gadget-2 and Pressure-Energy) are currently implemented.
+This must be directly implemented into the hydro schemes, therefore only a subset of schemes (Gadget-2, SPHENIX and Pressure-Energy) have the floor available.
 In order to implement it, you need equation 12 in `Hopkins 2013 <https://arxiv.org/abs/1206.5006>`_:
 
 .. math::

examples/SubgridTests/PressureFloor/README

+This example consists in a uniform box of particles.
+Using the cooling and gravity, the temperature of the gas will decrease
+slowly and the gas will start to collapse due to the numerical noise.
+If a pressure / energy floor is applied, the box should remain uniform.
+
+At the end, a plot containing the distribution of densities is generated.
\ No newline at end of file

examples/SubgridTests/PressureFloor/getGlass.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_16.hdf5

examples/SubgridTests/PressureFloor/makeIC.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2020 loic hausammann (loic.hausammann@epfl.ch)
+#
+# 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/>.
+#
+##############################################################################
+
+import h5py
+import numpy as np
+
+# Generates a SWIFT IC file with a constant density and pressure
+
+# Parameters
+periodic = 1  # 1 For periodic box
+boxSize = 200  # 200 kiloparsec
+rho = 200  # Density in code units
+T = 3000  # Initial Temperature
+gamma = 5./3.  # Gas adiabatic index
+fileName = "pressureFloor.hdf5"
+# ---------------------------------------------------
+
+# defines some constants
+# need to be changed in plotTemperature.py too
+h_frac = 0.76
+mu = 4. / (1. + 3. * h_frac)
+
+m_h_cgs = 1.67e-24
+k_b_cgs = 1.38e-16
+
+# defines units
+unit_length = 3.0857e21  # kpc
+unit_mass = 2.0e33  # solar mass
+unit_time = 3.0857e16  # ~ Gyr
+
+# Read id, position and h from glass
+glass = h5py.File("glassCube_16.hdf5", "r")
+ids = glass["/PartType0/ParticleIDs"][:]
+pos = glass["/PartType0/Coordinates"][:, :] * boxSize
+h = glass["/PartType0/SmoothingLength"][:] * boxSize
+
+# Compute basic properties
+numPart = np.size(pos) // 3
+mass = boxSize**3 * rho / numPart
+internalEnergy = k_b_cgs * T * mu / ((gamma - 1.) * m_h_cgs)
+internalEnergy *= (unit_time / unit_length)**2
+
+# File
+f = h5py.File(fileName, 'w')
+
+# Header
+grp = f.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+
+# Runtime parameters
+grp = f.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+# Units
+grp = f.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = unit_length
+grp.attrs["Unit mass in cgs (U_M)"] = unit_mass
+grp.attrs["Unit time in cgs (U_t)"] = unit_time
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+# Particle group
+grp = f.create_group("/PartType0")
+
+v = np.zeros((numPart, 3))
+ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+
+m = np.full((numPart, 1), mass)
+ds = grp.create_dataset('Masses', (numPart, 1), 'f')
+ds[()] = m
+
+h = np.reshape(h, (numPart, 1))
+ds = grp.create_dataset('SmoothingLength', (numPart, 1), 'f')
+ds[()] = h
+
+u = np.full((numPart, 1), internalEnergy)
+ds = grp.create_dataset('InternalEnergy', (numPart, 1), 'f')
+ds[()] = u
+
+ids = np.reshape(ids, (numPart, 1))
+ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
+ds[()] = ids
+
+ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = pos
+
+f.close()
+
+print("Initial condition generated")

examples/SubgridTests/PressureFloor/plotDensity.py

+#!/usr/bin/env python
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# 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/>.
+#
+##############################################################################
+
+import matplotlib.pyplot as plt
+from h5py import File
+import sys
+
+N = 100
+filename = "pressureFloor_%04i.hdf5" % int(sys.argv[-1])
+
+f = File(filename, "r")
+
+rho = f["PartType0"]["Densities"][:]
+
+plt.hist(rho, 100, log=True)
+plt.xlabel("Density [M$_\odot$ kpc$^{-3}$]")
+plt.ylabel("Number of particles")
+plt.savefig("density.png", dpi=200)

examples/SubgridTests/PressureFloor/pressureFloor.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     2.0e33     # Solar masses
+  UnitLength_in_cgs:   3.0857e21  # parsecs
+  UnitVelocity_in_cgs: 1.0e5      # megameters per second
+  UnitCurrent_in_cgs:  1          # Amperes
+  UnitTemp_in_cgs:     1          # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   8e1  # The end time of the simulation (in internal units).
+  dt_min:     1e-10  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e1  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            pressureFloor # Common part of the name of output files
+  time_first:          0.         # Time of the first output (in internal units)
+  delta_time:          1e0       # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e0 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.5      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature: 10.       # Kelvin
+  
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./pressureFloor.hdf5     # The file to read
+  periodic:   1
+  
+# Dimensionless pre-factor for the time-step condition
+LambdaCooling:
+  lambda_nH2_cgs:              1e-15 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3])
+  cooling_tstep_mult:          1e5        # Dimensionless pre-factor for the time-step condition
+
+# Parameters related to pressure floors    ----------------------------------------------
+GEARPressureFloor:
+  Jeans_factor: 10
+
+  # Parameters for the self-gravity scheme
+Gravity:
+  mesh_side_length:              12       # Number of cells along each axis for the periodic gravity mesh.
+  eta:                           2.5     # Constant dimensionless multiplier for time integration.
+  theta:                         0.7       # Opening angle (Multipole acceptance criterion).
+  max_physical_DM_softening:     0.5    # Maximal Plummer-equivalent softening length in physical coordinates for DM particles (in internal units).
+  max_physical_baryon_softening: 0.5    # 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.

examples/SubgridTests/PressureFloor/run.sh

+
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e glassCube_32.hdf5 ]
+then
+    echo "Fetching initial glass file for the cooling box example..."
+    ./getGlass.sh
+fi
+if [ ! -e coolingBox.hdf5 ]
+then
+    echo "Generating initial conditions for the cooling box example..."
+    python makeIC.py
+fi
+
+rm pressureFloor_*
+# Run SWIFT
+../../swift --self-gravity --hydro --cooling --threads=8 pressureFloor.yml
+
+# Check if the simulation collapsed
+python plotDensity.py 80

src/hydro/SPHENIX/hydro.h

@@ -37,6 +37,7 @@
 #include "hydro_space.h"
 #include "kernel_hydro.h"
 #include "minmax.h"
+#include "pressure_floor.h"
 
 #include "./hydro_parameters.h"
 
@@ -407,11 +408,14 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
 
   /* Compute the sound speed */
   const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
-  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+  const float pressure_including_floor =
+      pressure_floor_get_comoving_pressure(p, pressure, cosmo);
+  const float soundspeed =
+      gas_soundspeed_from_pressure(p->rho, pressure_including_floor);
 
   /* Update variables. */
   p->force.soundspeed = soundspeed;
-  p->force.pressure = pressure;
+  p->force.pressure = pressure_including_floor;
 
   p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
 }
@@ -589,9 +593,12 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
   /* Compute the norm of div v */
   const float abs_div_v = fabsf(p->viscosity.div_v);
 
-  /* Compute the sound speed -- see theory section for justification */
-  const float soundspeed = hydro_get_comoving_soundspeed(p);
+  /* Compute the sound speed  */
   const float pressure = hydro_get_comoving_pressure(p);
+  const float pressure_including_floor =
+      pressure_floor_get_comoving_pressure(p, pressure, cosmo);
+  const float soundspeed =
+      gas_soundspeed_from_pressure(p->rho, pressure_including_floor);
 
   /* Compute the Balsara switch */
   const float balsara =
@@ -610,7 +617,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
 
   /* Update variables. */
   p->force.f = grad_h_term;
-  p->force.pressure = pressure;
+  p->force.pressure = pressure_including_floor;
   p->force.soundspeed = soundspeed;
   p->force.balsara = balsara;
 }
@@ -719,7 +726,12 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float kernel_support_physical = p->h * cosmo->a * kernel_gamma;
   const float kernel_support_physical_inv = 1.f / kernel_support_physical;
   const float v_sig_physical = p->viscosity.v_sig * cosmo->a_factor_sound_speed;
-  const float soundspeed_physical = hydro_get_physical_soundspeed(p, cosmo);
+  const float pressure = hydro_get_comoving_pressure(p);
+  const float pressure_including_floor =
+      pressure_floor_get_comoving_pressure(p, pressure, cosmo);
+  const float soundspeed_physical =
+      gas_soundspeed_from_pressure(p->rho, pressure_including_floor) *
+      cosmo->a_factor_sound_speed;
 
   const float sound_crossing_time_inverse =
       soundspeed_physical * kernel_support_physical_inv;
@@ -850,9 +862,12 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
 
   /* Compute the sound speed */
   const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
-  const float soundspeed = hydro_get_comoving_soundspeed(p);
+  const float pressure_including_floor =
+      pressure_floor_get_comoving_pressure(p, pressure, cosmo);
+  const float soundspeed =
+      gas_soundspeed_from_pressure(p->rho, pressure_including_floor);
 
-  p->force.pressure = pressure;
+  p->force.pressure = pressure_including_floor;
   p->force.soundspeed = soundspeed;
 
   /* Update the signal velocity, if we need to. */
@@ -918,9 +933,12 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 
   /* Compute the new sound speed */
   const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
-  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+  const float pressure_including_floor =
+      pressure_floor_get_comoving_pressure(p, pressure, cosmo);
+  const float soundspeed =
+      gas_soundspeed_from_pressure(p->rho, pressure_including_floor);
 
-  p->force.pressure = pressure;
+  p->force.pressure = pressure_including_floor;
   p->force.soundspeed = soundspeed;
 
   /* Update signal velocity if we need to */
@@ -1033,9 +1051,12 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
   p->diffusion.alpha = hydro_props->diffusion.alpha;
 
   const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
-  const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
+  const float pressure_including_floor =
+      pressure_floor_get_comoving_pressure(p, pressure, cosmo);
+  const float soundspeed =
+      gas_soundspeed_from_pressure(p->rho, pressure_including_floor);
 
-  p->force.pressure = pressure;
+  p->force.pressure = pressure_including_floor;
   p->force.soundspeed = soundspeed;
 }
 

src/hydro/SPHENIX/hydro_part.h

@@ -29,6 +29,7 @@
 #include "black_holes_struct.h"
 #include "chemistry_struct.h"
 #include "cooling_struct.h"
+#include "pressure_floor_struct.h"
 #include "star_formation_struct.h"
 #include "timestep_limiter_struct.h"
 #include "tracers_struct.h"
@@ -197,6 +198,9 @@ struct part {
   /*! Black holes information (e.g. swallowing ID) */
   struct black_holes_part_data black_holes_data;
 
+  /* Additional data used by the pressure floor */
+  struct pressure_floor_part_data pressure_floor_data;
+
   /*! Time-step length */
   timebin_t time_bin;
 

src/pressure_floor.h

@@ -29,7 +29,6 @@
 
 /* Local includes */
 #include "common_io.h"
-#include "cosmology.h"
 #include "error.h"
 #include "inline.h"
 
@@ -37,7 +36,7 @@ extern struct pressure_floor_properties pressure_floor_props;
 
 /* Check if pressure floor is implemented in hydro */
 #ifndef PRESSURE_FLOOR_NONE
-#if defined(GADGET2_SPH) || defined(HOPKINS_PU_SPH)
+#if defined(GADGET2_SPH) || defined(HOPKINS_PU_SPH) || defined(SPHENIX_SPH)
 /* Implemented */
 #else
 #error Pressure floor not implemented with this hydro scheme

src/pressure_floor/GEAR/pressure_floor.h

@@ -20,6 +20,7 @@
 #define SWIFT_PRESSURE_FLOOR_GEAR_H
 
 /* Forward declaration */
+struct cosmology;
 __attribute__((always_inline)) static INLINE float
 pressure_floor_get_comoving_pressure(const struct part* p, const float pressure,
                                      const struct cosmology* cosmo);
@@ -160,8 +161,10 @@ __attribute__((always_inline)) INLINE static void pressure_floor_print_snapshot(
   io_write_attribute_s(h_grp, "Pressure floor", "GEAR");
 }
 
+#endif
+
 /**
- * @brief Finishes the density calculation.
+ * @brief Finishes the density calculation for the pressure floor properties.
  *
  * @param p The particle to act upon
  * @param cosmo The current cosmological model.
@@ -178,7 +181,8 @@ __attribute__((always_inline)) INLINE static void pressure_floor_end_density(
 }
 
 /**
- * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ * @brief Sets all the pressure floor fields to sensible values when the #part
+ * has 0 ngbs.
  *
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
@@ -225,5 +229,4 @@ pressure_floor_first_init_part(const struct phys_const* restrict phys_const,
   pressure_floor_init_part(p, xp);
 }
 
-#endif
 #endif /* SWIFT_PRESSURE_FLOOR_GEAR_H */

src/pressure_floor/none/pressure_floor.h

@@ -19,13 +19,14 @@
 #ifndef SWIFT_PRESSURE_FLOOR_NONE_H
 #define SWIFT_PRESSURE_FLOOR_NONE_H
 
-#include "adiabatic_index.h"
-#include "cosmology.h"
-#include "equation_of_state.h"
-#include "hydro_properties.h"
-#include "parser.h"
-#include "part.h"
-#include "units.h"
+/* Pre-declarations */
+struct cosmology;
+struct hydro_props;
+struct phys_const;
+struct part;
+struct xpart;
+struct swift_params;
+struct unit_system;
 
 /**
  * @file src/pressure_floor/none/pressure_floor.h
@@ -42,6 +43,9 @@ struct pressure_floor_properties {};
  *
  * Note that the particle is not updated!!
  *
+ * Since this is the 'none' model for floor, there is no floor and
+ * we just return the physical pressure that was received.
+ *
  * @param p The #part.
  * @param physical_pressure The physical pressure without any pressure floor.
  * @param cosmo The #cosmology model.
@@ -59,6 +63,9 @@ static INLINE float pressure_floor_get_physical_pressure(
  *
  * Note that the particle is not updated!!
  *
+ * Since this is the 'none' model for floor, there is no floor and
+ * we just return the comoving pressure that was received.
+ *
  * @param p The #part.
  * @param comoving_pressure The comoving pressure without any pressure floor.
  * @param cosmo The #cosmology model.
@@ -75,8 +82,7 @@ static INLINE float pressure_floor_get_comoving_pressure(
  * @brief Initialise the pressure floor by reading the parameters and converting
  * to internal units.
  *
- * The input temperatures and number densities are converted to pressure and
- * density assuming a neutral gas of primoridal abundance.
+ * Nothing to do here.
  *
  * @param params The YAML parameter file.
  * @param us The system of units used internally.
@@ -96,7 +102,10 @@ static INLINE void pressure_floor_init(struct pressure_floor_properties* props,
  * @param props The pressure floor properties.
  */
 static INLINE void pressure_floor_print(
-    const struct pressure_floor_properties* props) {}
+    const struct pressure_floor_properties* props) {
+
+  message("Pressure floor is 'none'");
+}
 
 #ifdef HAVE_HDF5
 
@@ -111,7 +120,9 @@ INLINE static void pressure_floor_print_snapshot(hid_t h_grp) {
 #endif
 
 /**
- * @brief Finishes the density calculation.
+ * @brief Finishes the density calculation for the pressure floor properties.
+ *
+ * Nothing to do here.
  *
  * @param p The particle to act upon
  * @param cosmo The current cosmological model.
@@ -120,7 +131,10 @@ __attribute__((always_inline)) INLINE static void pressure_floor_end_density(
     struct part* restrict p, const struct cosmology* cosmo) {}
 
 /**
- * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ * @brief Sets all the pressure floor fields to sensible values when the #part
+ * has 0 ngbs.
+ *
+ * Nothing to do here.
  *
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
@@ -135,6 +149,8 @@ pressure_floor_part_has_no_neighbours(struct part* restrict p,
  * @brief Sets the pressure_floor properties of the (x-)particles to a valid
  * start state.
  *
+ * Nothing to do here.
+ *
  * @param p Pointer to the particle data.
  * @param xp Pointer to the extended particle data.
  */
@@ -145,6 +161,8 @@ __attribute__((always_inline)) INLINE static void pressure_floor_init_part(
  * @brief Sets the pressure_floor properties of the (x-)particles to a valid
  * start state.
  *
+ * Nothing to do here.
+ *
  * @param phys_const The physical constant in internal units.
  * @param us The unit system.
  * @param cosmo The current cosmological model.