diff --git a/.gitignore b/.gitignore
index c70f1ab3aec9206e6c4465e736b529285daab5da..4a8a0b583aec90cffe3e450155fa54a4885778f0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -33,7 +33,12 @@ examples/*/*/*.xmf
examples/*/*/*.png
examples/*/*/*.txt
examples/*/*/*.dot
+examples/*/*/*.rst
+examples/*/*/*.hdf5
+examples/*/snapshots*
+examples/*/restart
examples/*/*/used_parameters.yml
+examples/*/*.mpg
examples/*/gravity_checks_*.dat
tests/testActivePair
@@ -59,7 +64,7 @@ tests/swift_dopair_125_standard.dat
tests/brute_force_125_perturbed.dat
tests/swift_dopair_125_perturbed.dat
tests/brute_force_pair_active.dat
-tests/brute_force_dopair2_active.dat
+tests/brute_force_dopair2_active.dat
tests/swift_dopair2_force_active.dat
tests/brute_force_periodic_BC_perturbed.dat
tests/swift_dopair_active.dat
@@ -108,6 +113,9 @@ tests/benchmarkInteractions
tests/testGravityDerivatives
tests/testPotentialSelf
tests/testPotentialPair
+tests/testEOS
+tests/testEOS*.txt
+tests/testEOS*.png
theory/latex/swift.pdf
theory/SPH/Kernels/kernels.pdf
diff --git a/configure.ac b/configure.ac
index dc0f3cb96526adb93a0291e5470113e5740b5df3..e3a796c14251e6a0309a04eab571a9e026c503e5 100644
--- a/configure.ac
+++ b/configure.ac
@@ -75,7 +75,7 @@ case "$with_subgrid" in
with_subgrid_chemistry=EAGLE
with_subgrid_hydro=gadget2
;;
- *)
+ *)
AC_MSG_ERROR([Unknown subgrid choice: $with_subgrid])
;;
esac
@@ -618,7 +618,7 @@ if test "x$with_fftw" != "xno"; then
AC_DEFINE([HAVE_FFTW],1,[The FFTW library appears to be present.]),
AC_MSG_ERROR(something is wrong with the FFTW library!), $FFTW_LIBS)
have_fftw="yes"
- fi
+ fi
if test "$have_fftw" = "no"; then
FFTW_LIBS=""
FFTW_INCS=""
@@ -855,7 +855,7 @@ fi
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
- [Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo, shadowfax debug default: gadget2@:>@]
+ [Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo, shadowfax, minimal-multi-mat, debug default: gadget2@:>@]
)],
[with_hydro="$withval"],
[with_hydro="gadget2"]
@@ -888,6 +888,9 @@ case "$with_hydro" in
shadowfax)
AC_DEFINE([SHADOWFAX_SPH], [1], [Shadowfax SPH])
;;
+ minimal-multi-mat)
+ AC_DEFINE([MINIMAL_MULTI_MAT_SPH], [1], [Minimal Multiple Material SPH])
+ ;;
*)
AC_MSG_ERROR([Unknown hydrodynamics scheme: $with_hydro])
@@ -975,7 +978,7 @@ esac
# Equation of state
AC_ARG_WITH([equation-of-state],
[AS_HELP_STRING([--with-equation-of-state=<EoS>],
- [equation of state @<:@ideal-gas, isothermal-gas default: ideal-gas@:>@]
+ [equation of state @<:@ideal-gas, isothermal-gas, planetary default: ideal-gas@:>@]
)],
[with_eos="$withval"],
[with_eos="ideal-gas"]
@@ -987,6 +990,9 @@ case "$with_eos" in
isothermal-gas)
AC_DEFINE([EOS_ISOTHERMAL_GAS], [1], [Isothermal gas equation of state])
;;
+ planetary)
+ AC_DEFINE([EOS_PLANETARY], [1], [All planetary equations of state])
+ ;;
*)
AC_MSG_ERROR([Unknown equation of state: $with_eos])
;;
@@ -1074,7 +1080,7 @@ case "$with_cooling" in
grackle)
AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
AC_DEFINE([COOLING_GRACKLE_MODE], [0], [Grackle chemistry network, mode 0])
- ;;
+ ;;
grackle1)
AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
AC_DEFINE([COOLING_GRACKLE_MODE], [1], [Grackle chemistry network, mode 1])
@@ -1082,11 +1088,11 @@ case "$with_cooling" in
grackle2)
AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
AC_DEFINE([COOLING_GRACKLE_MODE], [2], [Grackle chemistry network, mode 2])
- ;;
+ ;;
grackle3)
AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
AC_DEFINE([COOLING_GRACKLE_MODE], [3], [Grackle chemistry network, mode 3])
- ;;
+ ;;
EAGLE)
AC_DEFINE([COOLING_EAGLE], [1], [Cooling following the EAGLE model])
;;
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index 844a6fa68029af3c119d7c03e09e6dec957c5da4..a4a28480a09baf4ef0ef58a1bff3a537326e5f16 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -765,6 +765,7 @@ INPUT += @top_srcdir@/src/gravity/Default
INPUT += @top_srcdir@/src/stars/Default
INPUT += @top_srcdir@/src/riemann
INPUT += @top_srcdir@/src/potential/point_mass
+INPUT += @top_srcdir@/src/equation_of_state/ideal_gas
INPUT += @top_srcdir@/src/cooling/none
INPUT += @top_srcdir@/src/chemistry/none
diff --git a/examples/MoonFormingImpact/README.md b/examples/MoonFormingImpact/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..97a84f67c6aeeff4176a1385381f1cfe9e340c91
--- /dev/null
+++ b/examples/MoonFormingImpact/README.md
@@ -0,0 +1,34 @@
+Canonical Moon-Forming Giant Impact
+===================================
+
+NOTE: This doesn't really work because the EOS are different to Canup (2004) so
+the impactor just glances then flies away!
+
+A version of the canonical moon-forming giant impact of Theia onto the early
+Earth (Canup 2004; Barr 2016). Both bodies are here made of a (Tillotson) iron
+core and granite mantle. Only ~10,000 particles are used for a quick and crude
+simulation.
+
+Setup
+-----
+
+In `swiftsim/`:
+
+`$ ./configure --with-hydro=minimal-multi-mat --with-equation-of-state=planetary`
+`$ make`
+
+In `swiftsim/examples/MoonFormingImpact/`:
+
+`$ ./get_init_cond.sh`
+
+Run
+---
+
+`$ ./run.sh`
+
+Output
+------
+
+`$ python plot.py`
+`$ mplayer anim.mpg`
+
diff --git a/examples/MoonFormingImpact/get_init_cond.sh b/examples/MoonFormingImpact/get_init_cond.sh
new file mode 100755
index 0000000000000000000000000000000000000000..7d63943c2c5dc3bd4ab88e63a2abba62cc3f04a5
--- /dev/null
+++ b/examples/MoonFormingImpact/get_init_cond.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/moon_forming_impact.hdf5
diff --git a/examples/MoonFormingImpact/moon_forming_impact.yml b/examples/MoonFormingImpact/moon_forming_impact.yml
new file mode 100644
index 0000000000000000000000000000000000000000..7d726bc02ba4fb6c8dd02f74907ffc48c3ed9431
--- /dev/null
+++ b/examples/MoonFormingImpact/moon_forming_impact.yml
@@ -0,0 +1,44 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 5.9724e27 # Grams
+ UnitLength_in_cgs: 6.371e8 # Centimeters
+ UnitVelocity_in_cgs: 6.371e8 # Centimeters 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: 100000 # The end time of the simulation (in internal units).
+ dt_min: 0.001 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 100 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ # Common part of the name of output files
+ basename: snapshots/moon_forming_impact
+ time_first: 0 # Time of the first output (in internal units)
+ delta_time: 100 # Time difference between consecutive outputs (in internal units)
+ label_delta: 100 # Integer increment between snapshot output labels
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 500 # 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).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.2 # Courant-Friedrich-Levy condition for time integration.
+
+# 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.005 # Comoving softening length (in internal units).
+ max_physical_softening: 0.005 # Physical softening length (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+ # The initial conditions file to read
+ file_name: moon_forming_impact.hdf5
diff --git a/examples/MoonFormingImpact/plot.py b/examples/MoonFormingImpact/plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..aa0d64a5d0d06709d51b1db231c507e22861f36c
--- /dev/null
+++ b/examples/MoonFormingImpact/plot.py
@@ -0,0 +1,285 @@
+"""
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 Jacob Kegerreis (jacob.kegerreis@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/>.
+#
+###############################################################################
+
+Plotting script for the Canonical Moon-Forming Giant Impact example.
+
+Save a figure for each snapshot in `./plots/` then make them into a simple
+animation with ffmpeg in `./`.
+
+Usage:
+ `$ python plot.py time_end delta_time`
+
+ Sys args:
+ + `time_end` | (opt) int | The time of the last snapshot to plot.
+ Default = 100000
+ + `delta_time` | (opt) int | The time between successive snapshots.
+ Default = 100
+"""
+
+from __future__ import division
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+import h5py
+import sys
+import subprocess
+
+# Particle array fields
+dtype_picle = [
+ ('m', float), ('x', float), ('y', float), ('z', float), ('v_x', float),
+ ('v_y', float), ('v_z', float), ('ID', int), ('rho', float), ('u', float),
+ ('phi', float), ('P', float), ('h', float), ('mat_ID', int), ('r', float)
+ ]
+
+s_to_hour = 1 / 60**2
+
+# Snapshot info
+file_snap = "./snapshots/moon_forming_impact_"
+file_plot = "./plots/moon_forming_impact_"
+# Number of particles in the target body
+num_target = 9496
+
+# Material types (copied from src/equation_of_state/planetary/equation_of_state.h)
+type_factor = 100
+Di_type = {
+ 'Til' : 1,
+ 'HM80' : 2,
+ 'ANEOS' : 3,
+ 'SESAME' : 4,
+}
+Di_material = {
+ # Tillotson
+ 'Til_iron' : Di_type['Til']*type_factor,
+ 'Til_granite' : Di_type['Til']*type_factor + 1,
+ 'Til_water' : Di_type['Til']*type_factor + 2,
+ # Hubbard & MacFarlane (1980) Uranus/Neptune
+ 'HM80_HHe' : Di_type['HM80']*type_factor, # Hydrogen-helium atmosphere
+ 'HM80_ice' : Di_type['HM80']*type_factor + 1, # H20-CH4-NH3 ice mix
+ 'HM80_rock' : Di_type['HM80']*type_factor + 2, # SiO2-MgO-FeS-FeO rock mix
+ # ANEOS
+ 'ANEOS_iron' : Di_type['ANEOS']*type_factor,
+ 'MANEOS_forsterite' : Di_type['ANEOS']*type_factor + 1,
+ # SESAME
+ 'SESAME_iron' : Di_type['SESAME']*type_factor,
+}
+
+# Material offset for impactor particles
+ID_imp = 10000
+# Material colours
+Di_mat_colour = {
+ # Target
+ Di_material['Til_iron'] : 'orange',
+ Di_material['Til_granite'] : '#FFF0E0',
+ # Impactor
+ Di_material['Til_iron'] + ID_imp : 'dodgerblue',
+ Di_material['Til_granite'] + ID_imp : '#A080D0',
+ }
+
+
+def load_snapshot(filename):
+ """ Load the hdf5 snapshot file and return the structured particle array.
+ """
+ # Add extension if needed
+ if (filename[-5:] != ".hdf5"):
+ filename += ".hdf5"
+
+ # Load the hdf5 file
+ with h5py.File(filename, 'r') as f:
+ header = f['Header'].attrs
+ A2_pos = f['PartType0/Coordinates'].value
+ A2_vel = f['PartType0/Velocities'].value
+
+ # Structured array of all particle data
+ A2_picle = np.empty(header['NumPart_Total'][0],
+ dtype=dtype_picle)
+
+ A2_picle['x'] = A2_pos[:, 0]
+ A2_picle['y'] = A2_pos[:, 1]
+ A2_picle['z'] = A2_pos[:, 2]
+ A2_picle['v_x'] = A2_vel[:, 0]
+ A2_picle['v_y'] = A2_vel[:, 1]
+ A2_picle['v_z'] = A2_vel[:, 2]
+ A2_picle['m'] = f['PartType0/Masses'].value
+ A2_picle['ID'] = f['PartType0/ParticleIDs'].value
+ A2_picle['rho'] = f['PartType0/Density'].value
+ A2_picle['u'] = f['PartType0/InternalEnergy'].value
+ A2_picle['phi'] = f['PartType0/Potential'].value
+ A2_picle['P'] = f['PartType0/Pressure'].value
+ A2_picle['h'] = f['PartType0/SmoothingLength'].value
+ A2_picle['mat_ID'] = f['PartType0/MaterialID'].value
+
+ return A2_picle
+
+
+def process_particles(A2_picle, num_target):
+ """ Modify things like particle units, material IDs, and coordinate origins.
+ """
+ # Offset material IDs for impactor particles
+ A2_picle['mat_ID'][A2_picle['ID'] >= num_target] += ID_imp
+
+ # Shift coordinates to the centre of the target's core's mass and momentum
+ sel_tar = np.where(A2_picle['mat_ID'] == Di_material['Til_iron'])[0]
+
+ # Centre of mass
+ m_tot = np.sum(A2_picle[sel_tar]['m'])
+ x_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['x']) / m_tot
+ y_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['y']) / m_tot
+ z_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['z']) / m_tot
+
+ # Change origin to the centre-of-mass
+ A2_picle['x'] -= x_com
+ A2_picle['y'] -= y_com
+ A2_picle['z'] -= z_com
+ A2_picle['r'] = np.sqrt(
+ A2_picle['x']**2 + A2_picle['y']**2 + A2_picle['z']**2
+ )
+
+ # Centre of momentum
+ v_x_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['v_x']) / m_tot
+ v_y_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['v_y']) / m_tot
+ v_z_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['v_z']) / m_tot
+
+ # Change to the centre-of-momentum frame of reference
+ A2_picle['v_x'] -= v_x_com
+ A2_picle['v_y'] -= v_y_com
+ A2_picle['v_z'] -= v_z_com
+
+ return A2_picle
+
+
+def plot_snapshot(A2_picle, filename, time, ax_lim=100, dz=0.1):
+ """ Plot the snapshot particles and save the figure.
+ """
+ # Add extension if needed
+ if (filename[-5:] != ".png"):
+ filename += ".png"
+
+ fig = plt.figure(figsize=(9, 9))
+ ax = fig.add_subplot(111, aspect='equal')
+
+ # Plot slices in z below zero
+ for z in np.arange(-ax_lim, 0, dz):
+ sel_z = np.where((z < A2_picle['z']) & (A2_picle['z'] < z+dz))[0]
+ A2_picle_z = A2_picle[sel_z]
+
+ # Plot each material
+ for mat_ID, colour in Di_mat_colour.iteritems():
+ sel_col = np.where(A2_picle_z['mat_ID'] == mat_ID)[0]
+
+ ax.scatter(
+ A2_picle_z[sel_col]['x'], A2_picle_z[sel_col]['y'],
+ c=colour, edgecolors='none', marker='.', s=50, alpha=0.7
+ )
+
+ # Axes etc.
+ ax.set_axis_bgcolor('k')
+
+ ax.set_xlabel("x Position ($R_\oplus$)")
+ ax.set_ylabel("y Position ($R_\oplus$)")
+
+ ax.set_xlim(-ax_lim, ax_lim)
+ ax.set_ylim(-ax_lim, ax_lim)
+
+ plt.text(
+ -0.92*ax_lim, 0.85*ax_lim, "%.1f h" % (time*s_to_hour), fontsize=20,
+ color='w'
+ )
+
+ # Font sizes
+ for item in (
+ [ax.title, ax.xaxis.label, ax.yaxis.label] + ax.get_xticklabels() +
+ ax.get_yticklabels()
+ ):
+ item.set_fontsize(20)
+
+ plt.tight_layout()
+
+ plt.savefig(filename)
+ plt.close()
+
+
+if __name__ == '__main__':
+ # Sys args
+ try:
+ time_end = int(sys.argv[1])
+
+ try:
+ delta_time = int(sys.argv[2])
+ except IndexError:
+ delta_time = 100
+ except IndexError:
+ time_end = 100000
+ delta_time = 100
+
+ # Load and plot each snapshot
+ for i_snap in range(int(time_end/delta_time) + 1):
+ snap_time = i_snap * delta_time
+ print "\rPlotting snapshot %06d (%d of %d)" % (
+ snap_time, i_snap+1, int(time_end/delta_time)
+ ),
+ sys.stdout.flush()
+
+ # Load particle data
+ filename = "%s%06d" % (file_snap, snap_time)
+ A2_picle = load_snapshot(filename)
+
+ # Process particle data
+ A2_picle = process_particles(A2_picle, num_target)
+
+ # Plot particles
+ filename = "%s%06d" % (file_plot, snap_time)
+ plot_snapshot(A2_picle, filename, snap_time)
+
+ # Animation
+ command = (
+ "ffmpeg -framerate 12 -i plots/moon_forming_impact_%*.png -r 25 "
+ "anim.mpg -y"
+ )
+ print "\n%s\n" % command
+ subprocess.check_output(command, shell=True)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/examples/MoonFormingImpact/run.sh b/examples/MoonFormingImpact/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..165dae3a24a9c30960959fbb37aa6e1da2eb851f
--- /dev/null
+++ b/examples/MoonFormingImpact/run.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+../swift -G -s -t 8 moon_forming_impact.yml
diff --git a/examples/UranusImpact/README.md b/examples/UranusImpact/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..178a3937ecbe527df8e8e82a0d8fd8bcbf9dbef7
--- /dev/null
+++ b/examples/UranusImpact/README.md
@@ -0,0 +1,40 @@
+Uranus Giant Impact
+===================
+
+A simple version of the low angular momentum impact onto the early Uranus shown
+in Kegerreis et al. (2018), Fig. 2; with only ~10,000 particles for a quick and
+crude simulation.
+
+The collision of a 2 Earth mass impactor onto a proto-Uranus that can explain
+the spin of the present-day planet, with an angular momentum of 2e36 kg m^2 s^-1
+and velocity at infinity of 5 km s^-1 for a relatively head-on impact.
+
+Both bodies have a rocky core and icy mantle, with a hydrogen-helium atmosphere
+on the target as well. Although with this low number of particles it cannot be
+modelled in any detail.
+
+Setup
+-----
+
+In `swiftsim/`:
+
+`$ ./configure --with-hydro=minimal-multi-mat --with-equation-of-state=planetary`
+
+`$ make`
+
+In `swiftsim/examples/UranusImpact/`:
+
+`$ ./get_init_cond.sh`
+
+Run
+---
+
+`$ ./run.sh`
+
+Analysis
+--------
+
+`$ python plot.py`
+
+`$ mplayer anim.mpg`
+
diff --git a/examples/UranusImpact/get_init_cond.sh b/examples/UranusImpact/get_init_cond.sh
new file mode 100755
index 0000000000000000000000000000000000000000..e12e009adfbd727cb2452ac21c477b3ecd77b9c9
--- /dev/null
+++ b/examples/UranusImpact/get_init_cond.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/uranus_impact.hdf5
diff --git a/examples/UranusImpact/plot.py b/examples/UranusImpact/plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..3db3bf21bb15862ec524a069c38e47564b48df1d
--- /dev/null
+++ b/examples/UranusImpact/plot.py
@@ -0,0 +1,291 @@
+"""
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 Jacob Kegerreis (jacob.kegerreis@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/>.
+#
+###############################################################################
+
+Plotting script for the Uranus Giant Impact example.
+
+Save a figure for each snapshot in `./plots/` then make them into a simple
+animation with ffmpeg in `./`.
+
+The snapshot plots show all particles with z < 0, coloured by their material.
+
+Usage:
+ `$ python plot.py time_end delta_time`
+
+ Sys args:
+ + `time_end` | (opt) int | The time of the last snapshot to plot.
+ Default = 100000
+ + `delta_time` | (opt) int | The time between successive snapshots.
+ Default = 500
+"""
+
+from __future__ import division
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+import h5py
+import sys
+import subprocess
+
+# Particle array fields
+dtype_picle = [
+ ('m', float), ('x', float), ('y', float), ('z', float), ('v_x', float),
+ ('v_y', float), ('v_z', float), ('ID', int), ('rho', float), ('u', float),
+ ('phi', float), ('P', float), ('h', float), ('mat_ID', int), ('r', float)
+ ]
+
+s_to_hour = 1 / 60**2
+R_Ea = 6.371e6
+
+# Default sys args
+time_end_default = 100000
+delta_time_default = 500
+
+# Snapshot info
+file_snap = "./snapshots/uranus_impact_"
+file_plot = "./plots/uranus_impact_"
+
+# Number of particles in the target body
+num_target = 8992
+
+# Material types (copied from src/equation_of_state/planetary/equation_of_state.h)
+type_factor = 100
+Di_type = {
+ 'Til' : 1,
+ 'HM80' : 2,
+ 'ANEOS' : 3,
+ 'SESAME' : 4,
+}
+Di_material = {
+ # Tillotson
+ 'Til_iron' : Di_type['Til']*type_factor,
+ 'Til_granite' : Di_type['Til']*type_factor + 1,
+ 'Til_water' : Di_type['Til']*type_factor + 2,
+ # Hubbard & MacFarlane (1980) Uranus/Neptune
+ 'HM80_HHe' : Di_type['HM80']*type_factor, # Hydrogen-helium atmosphere
+ 'HM80_ice' : Di_type['HM80']*type_factor + 1, # H20-CH4-NH3 ice mix
+ 'HM80_rock' : Di_type['HM80']*type_factor + 2, # SiO2-MgO-FeS-FeO rock mix
+ # ANEOS
+ 'ANEOS_iron' : Di_type['ANEOS']*type_factor,
+ 'MANEOS_forsterite' : Di_type['ANEOS']*type_factor + 1,
+ # SESAME
+ 'SESAME_iron' : Di_type['SESAME']*type_factor,
+}
+
+# Material offset for impactor particles
+ID_imp = 10000
+# Material colours
+Di_mat_colour = {
+ # Target
+ Di_material['HM80_HHe'] : '#33DDFF',
+ Di_material['HM80_ice'] : 'lightsteelblue',
+ Di_material['HM80_rock'] : 'slategrey',
+ # Impactor
+ Di_material['HM80_ice'] + ID_imp : '#A080D0',
+ Di_material['HM80_rock'] + ID_imp : '#706050',
+ }
+
+
+def load_snapshot(filename):
+ """ Load the hdf5 snapshot file and return the structured particle array.
+ """
+ # Add extension if needed
+ if (filename[-5:] != ".hdf5"):
+ filename += ".hdf5"
+
+ # Load the hdf5 file
+ with h5py.File(filename, 'r') as f:
+ header = f['Header'].attrs
+ A2_pos = f['PartType0/Coordinates'].value
+ A2_vel = f['PartType0/Velocities'].value
+
+ # Structured array of all particle data
+ A2_picle = np.empty(header['NumPart_Total'][0],
+ dtype=dtype_picle)
+
+ A2_picle['x'] = A2_pos[:, 0]
+ A2_picle['y'] = A2_pos[:, 1]
+ A2_picle['z'] = A2_pos[:, 2]
+ A2_picle['v_x'] = A2_vel[:, 0]
+ A2_picle['v_y'] = A2_vel[:, 1]
+ A2_picle['v_z'] = A2_vel[:, 2]
+ A2_picle['m'] = f['PartType0/Masses'].value
+ A2_picle['ID'] = f['PartType0/ParticleIDs'].value
+ A2_picle['rho'] = f['PartType0/Density'].value
+ A2_picle['u'] = f['PartType0/InternalEnergy'].value
+ A2_picle['phi'] = f['PartType0/Potential'].value
+ A2_picle['P'] = f['PartType0/Pressure'].value
+ A2_picle['h'] = f['PartType0/SmoothingLength'].value
+ A2_picle['mat_ID'] = f['PartType0/MaterialID'].value
+
+ return A2_picle
+
+
+def process_particles(A2_picle, num_target):
+ """ Modify things like particle units, material IDs, and coordinate origins.
+ """
+ # Offset material IDs for impactor particles
+ A2_picle['mat_ID'][A2_picle['ID'] >= num_target] += ID_imp
+
+ # Shift coordinates to the centre of the target's core's mass and momentum
+ sel_tar = np.where(A2_picle['mat_ID'] == Di_material['HM80_rock'])[0]
+
+ # Centre of mass
+ m_tot = np.sum(A2_picle[sel_tar]['m'])
+ x_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['x']) / m_tot
+ y_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['y']) / m_tot
+ z_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['z']) / m_tot
+
+ # Change origin to the centre-of-mass
+ A2_picle['x'] -= x_com
+ A2_picle['y'] -= y_com
+ A2_picle['z'] -= z_com
+ A2_picle['r'] = np.sqrt(
+ A2_picle['x']**2 + A2_picle['y']**2 + A2_picle['z']**2
+ )
+
+ # Centre of momentum
+ v_x_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['v_x']) / m_tot
+ v_y_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['v_y']) / m_tot
+ v_z_com = np.sum(A2_picle[sel_tar]['m'] * A2_picle[sel_tar]['v_z']) / m_tot
+
+ # Change to the centre-of-momentum frame of reference
+ A2_picle['v_x'] -= v_x_com
+ A2_picle['v_y'] -= v_y_com
+ A2_picle['v_z'] -= v_z_com
+
+ return A2_picle
+
+
+def plot_snapshot(A2_picle, filename, time, ax_lim=13, dz=0.1):
+ """ Plot the snapshot particles and save the figure.
+ """
+ # Add extension if needed
+ if (filename[-5:] != ".png"):
+ filename += ".png"
+
+ fig = plt.figure(figsize=(9, 9))
+ ax = fig.add_subplot(111, aspect='equal')
+
+ # Plot slices in z below zero
+ for z in np.arange(-ax_lim, 0, dz):
+ sel_z = np.where((z < A2_picle['z']) & (A2_picle['z'] < z+dz))[0]
+ A2_picle_z = A2_picle[sel_z]
+
+ # Plot each material
+ for mat_ID, colour in Di_mat_colour.iteritems():
+ sel_col = np.where(A2_picle_z['mat_ID'] == mat_ID)[0]
+
+ ax.scatter(
+ A2_picle_z[sel_col]['x'], A2_picle_z[sel_col]['y'],
+ c=colour, edgecolors='none', marker='.', s=50, alpha=0.7
+ )
+
+ # Axes etc.
+ ax.set_axis_bgcolor('k')
+
+ ax.set_xlabel("x Position ($R_\oplus$)")
+ ax.set_ylabel("y Position ($R_\oplus$)")
+
+ ax.set_xlim(-ax_lim, ax_lim)
+ ax.set_ylim(-ax_lim, ax_lim)
+
+ plt.text(
+ -0.92*ax_lim, 0.85*ax_lim, "%.1f h" % (time*s_to_hour), fontsize=20,
+ color='w'
+ )
+
+ # Font sizes
+ for item in (
+ [ax.title, ax.xaxis.label, ax.yaxis.label] + ax.get_xticklabels() +
+ ax.get_yticklabels()
+ ):
+ item.set_fontsize(20)
+
+ plt.tight_layout()
+
+ plt.savefig(filename)
+ plt.close()
+
+
+if __name__ == '__main__':
+ # Sys args
+ try:
+ time_end = int(sys.argv[1])
+ try:
+ delta_time = int(sys.argv[2])
+ except IndexError:
+ delta_time = delta_time_default
+ except IndexError:
+ time_end = time_end_default
+ delta_time = delta_time_default
+
+ # Load and plot each snapshot
+ for i_snap in range(int(time_end/delta_time) + 1):
+ snap_time = i_snap * delta_time
+ print "\rPlotting snapshot %06d (%d of %d)" % (
+ snap_time, i_snap+1, int(time_end/delta_time)
+ ),
+ sys.stdout.flush()
+
+ # Load particle data
+ filename = "%s%06d" % (file_snap, snap_time)
+ A2_picle = load_snapshot(filename)
+
+ # Process particle data
+ A2_picle = process_particles(A2_picle, num_target)
+
+ # Plot particles
+ filename = "%s%06d" % (file_plot, snap_time)
+ plot_snapshot(A2_picle, filename, snap_time)
+
+ # Animation
+ command = (
+ "ffmpeg -framerate 10 -i plots/uranus_impact_%*.png -r 25 anim.mpg -y"
+ )
+ print "\n$ %s\n" % command
+ subprocess.call(command, shell=True)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/examples/UranusImpact/run.sh b/examples/UranusImpact/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..c6773b7e40fff3fa312dfcb5ba4ada9d9e4b1b8d
--- /dev/null
+++ b/examples/UranusImpact/run.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+../swift -G -s -t 8 uranus_impact.yml
diff --git a/examples/UranusImpact/uranus_impact.yml b/examples/UranusImpact/uranus_impact.yml
new file mode 100644
index 0000000000000000000000000000000000000000..aae9f66847a9f9b55984ea5d2ea1c79099c01e95
--- /dev/null
+++ b/examples/UranusImpact/uranus_impact.yml
@@ -0,0 +1,43 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 5.9724e27 # Grams
+ UnitLength_in_cgs: 6.371e8 # Centimeters
+ UnitVelocity_in_cgs: 6.371e8 # Centimeters 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: 100000 # The end time of the simulation (in internal units).
+ dt_min: 0.001 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 100 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ # Common part of the name of output files
+ basename: snapshots/uranus_impact
+ time_first: 0 # Time of the first output (in internal units)
+ delta_time: 500 # Time difference between consecutive outputs (in internal units)
+ label_delta: 500 # Integer increment between snapshot output labels
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1000 # 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).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.2 # Courant-Friedrich-Levy condition for time integration.
+
+# 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.01 # Comoving softening length (in internal units).
+ max_physical_softening: 0.01 # Physical softening length (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: uranus_impact.hdf5 # The initial conditions file to read
diff --git a/src/debug.c b/src/debug.c
index 4000625a5e7d1bc83600a3f20a0db160057f07a6..bf167e9a2954ed7b81588bc2aacf8e3b2b8cbe09 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -52,6 +52,8 @@
#include "./hydro/Gizmo/hydro_debug.h"
#elif defined(SHADOWFAX_SPH)
#include "./hydro/Shadowswift/hydro_debug.h"
+#elif defined(MINIMAL_MULTI_MAT_SPH)
+#include "./hydro/MinimalMultiMat/hydro_debug.h"
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/equation_of_state.h b/src/equation_of_state.h
index 195b52514f2acc0c40959e09c088a06f0a411869..1fed652f354e01a89f5e0b0540d9d326a7750b02 100644
--- a/src/equation_of_state.h
+++ b/src/equation_of_state.h
@@ -34,6 +34,8 @@
#include "./equation_of_state/ideal_gas/equation_of_state.h"
#elif defined(EOS_ISOTHERMAL_GAS)
#include "./equation_of_state/isothermal/equation_of_state.h"
+#elif defined(EOS_PLANETARY)
+#include "./equation_of_state/planetary/equation_of_state.h"
#else
#error "Invalid choice of equation of state"
#endif
diff --git a/src/equation_of_state/planetary/aneos.h b/src/equation_of_state/planetary/aneos.h
new file mode 100644
index 0000000000000000000000000000000000000000..904288b2fdf3ba825cdc7d114ebb61cd42de198d
--- /dev/null
+++ b/src/equation_of_state/planetary/aneos.h
@@ -0,0 +1,144 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * 2018 Jacob Kegerreis (jacob.kegerreis@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_ANEOS_EQUATION_OF_STATE_H
+#define SWIFT_ANEOS_EQUATION_OF_STATE_H
+
+/**
+ * @file equation_of_state/planetary/aneos.h
+ *
+ * Contains the (M)ANEOS EOS functions for
+ * equation_of_state/planetary/equation_of_state.h
+ *
+ * Adapted from the implementation in Gadget 2 of Cuk & Stewart (2012)
+ *
+ */
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "common_io.h"
+#include "equation_of_state.h"
+#include "inline.h"
+#include "physical_constants.h"
+#include "units.h"
+
+// ANEOS parameters
+struct ANEOS_params {
+ enum eos_planetary_material_id mat_id;
+};
+
+// Parameter values for each material (cgs units)
+INLINE static void set_ANEOS_iron(struct ANEOS_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ mat->mat_id = mat_id;
+}
+INLINE static void set_MANEOS_forsterite(
+ struct ANEOS_params *mat, enum eos_planetary_material_id mat_id) {
+ mat->mat_id = mat_id;
+}
+
+// Convert from cgs to internal units
+INLINE static void convert_units_ANEOS(struct ANEOS_params *mat,
+ const struct unit_system *us) {}
+
+// gas_internal_energy_from_entropy
+INLINE static float ANEOS_internal_energy_from_entropy(
+ float density, float entropy, const struct ANEOS_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_pressure_from_entropy
+INLINE static float ANEOS_pressure_from_entropy(
+ float density, float entropy, const struct ANEOS_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_entropy_from_pressure
+INLINE static float ANEOS_entropy_from_pressure(
+ float density, float pressure, const struct ANEOS_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_soundspeed_from_entropy
+INLINE static float ANEOS_soundspeed_from_entropy(
+ float density, float entropy, const struct ANEOS_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_entropy_from_internal_energy
+INLINE static float ANEOS_entropy_from_internal_energy(
+ float density, float u, const struct ANEOS_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_pressure_from_internal_energy
+INLINE static float ANEOS_pressure_from_internal_energy(
+ float density, float u, const struct ANEOS_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_internal_energy_from_pressure
+INLINE static float ANEOS_internal_energy_from_pressure(
+ float density, float P, const struct ANEOS_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_soundspeed_from_internal_energy
+INLINE static float ANEOS_soundspeed_from_internal_energy(
+ float density, float u, const struct ANEOS_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_soundspeed_from_pressure
+INLINE static float ANEOS_soundspeed_from_pressure(
+ float density, float P, const struct ANEOS_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+#endif /* SWIFT_ANEOS_EQUATION_OF_STATE_H */
diff --git a/src/equation_of_state/planetary/equation_of_state.h b/src/equation_of_state/planetary/equation_of_state.h
new file mode 100644
index 0000000000000000000000000000000000000000..fcc1c959152801253d465915d8b28b1552c913b8
--- /dev/null
+++ b/src/equation_of_state/planetary/equation_of_state.h
@@ -0,0 +1,1147 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * 2018 Jacob Kegerreis (jacob.kegerreis@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_PLANETARY_EQUATION_OF_STATE_H
+#define SWIFT_PLANETARY_EQUATION_OF_STATE_H
+
+/**
+ * @file equation_of_state/planetary/equation_of_state.h
+ *
+ * For any/all of the planetary EOS. Each EOS type's functions are set in its
+ * own header file: `equation_of_state/planetary/<eos_type>.h`.
+ * See `material_id` for the available choices.
+ *
+ * Not all functions are implemented for all EOS types, so not all can be used
+ * with all hydro formulations yet.
+ */
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "common_io.h"
+#include "inline.h"
+#include "physical_constants.h"
+#include "units.h"
+
+extern struct eos_parameters eos;
+
+/*! Material identifier flags (material_ID = type_ID * type_factor + unit_ID) */
+#define eos_planetary_type_factor 100
+
+/**
+ * @brief Master type for the planetary equation of state.
+ */
+enum eos_planetary_type_id {
+ eos_planetary_type_Til = 1,
+ eos_planetary_type_HM80 = 2,
+ eos_planetary_type_ANEOS = 3,
+ eos_planetary_type_SESAME = 4,
+} __attribute__((packed));
+
+/**
+ * @brief Minor type for the planetary equation of state.
+ */
+enum eos_planetary_material_id {
+
+ /* Tillotson */
+
+ /*! Tillotson iron */
+ eos_planetary_id_Til_iron =
+ eos_planetary_type_Til * eos_planetary_type_factor,
+
+ /*! Tillotson granite */
+ eos_planetary_id_Til_granite =
+ eos_planetary_type_Til * eos_planetary_type_factor + 1,
+
+ /*! Tillotson water */
+ eos_planetary_id_Til_water =
+ eos_planetary_type_Til * eos_planetary_type_factor + 2,
+
+ /* Hubbard & MacFarlane (1980) Uranus/Neptune */
+
+ /*! Hydrogen-helium atmosphere */
+ eos_planetary_id_HM80_HHe =
+ eos_planetary_type_HM80 * eos_planetary_type_factor,
+
+ /*! H20-CH4-NH3 ice mix */
+ eos_planetary_id_HM80_ice =
+ eos_planetary_type_HM80 * eos_planetary_type_factor + 1,
+
+ /*! SiO2-MgO-FeS-FeO rock mix */
+ eos_planetary_id_HM80_rock =
+ eos_planetary_type_HM80 * eos_planetary_type_factor + 2,
+
+ /* ANEOS */
+
+ /*! ANEOS iron */
+ eos_planetary_id_ANEOS_iron =
+ eos_planetary_type_ANEOS * eos_planetary_type_factor,
+
+ /*! ANEOS forsterite */
+ eos_planetary_id_MANEOS_forsterite =
+ eos_planetary_type_ANEOS * eos_planetary_type_factor + 1,
+
+ /* SESAME */
+
+ /*! SESAME iron */
+ eos_planetary_id_SESAME_iron =
+ eos_planetary_type_SESAME * eos_planetary_type_factor,
+} __attribute__((packed));
+
+/* Individual EOS function headers. */
+#include "aneos.h"
+#include "hm80.h"
+#include "sesame.h"
+#include "tillotson.h"
+
+/**
+ * @brief The parameters of the equation of state.
+ */
+struct eos_parameters {
+ struct Til_params Til_iron, Til_granite, Til_water;
+ struct HM80_params HM80_HHe, HM80_ice, HM80_rock;
+ struct ANEOS_params ANEOS_iron, MANEOS_forsterite;
+ struct SESAME_params SESAME_iron;
+};
+
+/**
+ * @brief Returns the internal energy given density and entropy
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_entropy(float density, float entropy,
+ enum eos_planetary_material_id mat_id) {
+
+ const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+
+ /* Select the material base type */
+ switch (type) {
+
+ /* Tillotson EoS */
+ case eos_planetary_type_Til:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ return Til_internal_energy_from_entropy(density, entropy,
+ &eos.Til_iron);
+ break;
+
+ case eos_planetary_id_Til_granite:
+ return Til_internal_energy_from_entropy(density, entropy,
+ &eos.Til_granite);
+ break;
+
+ case eos_planetary_id_Til_water:
+ return Til_internal_energy_from_entropy(density, entropy,
+ &eos.Til_water);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* Hubbard & MacFarlane (1980) EoS */
+ case eos_planetary_type_HM80:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ return HM80_internal_energy_from_entropy(density, entropy,
+ &eos.HM80_HHe);
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ return HM80_internal_energy_from_entropy(density, entropy,
+ &eos.HM80_ice);
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ return HM80_internal_energy_from_entropy(density, entropy,
+ &eos.HM80_rock);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* ANEOS EoS */
+ case eos_planetary_type_ANEOS:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_iron:
+ return ANEOS_internal_energy_from_entropy(density, entropy,
+ &eos.ANEOS_iron);
+ break;
+
+ case eos_planetary_id_MANEOS_forsterite:
+ return ANEOS_internal_energy_from_entropy(density, entropy,
+ &eos.MANEOS_forsterite);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_internal_energy_from_entropy(density, entropy,
+ &eos.SESAME_iron);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ default:
+ error("Unknown material type! mat_id = %d", mat_id);
+ return 0.f;
+ }
+}
+
+/**
+ * @brief Returns the pressure given density and entropy
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
+ float density, float entropy, enum eos_planetary_material_id mat_id) {
+
+ const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+
+ /* Select the material base type */
+ switch (type) {
+
+ /* Tillotson EoS */
+ case eos_planetary_type_Til:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ return Til_pressure_from_entropy(density, entropy, &eos.Til_iron);
+ break;
+
+ case eos_planetary_id_Til_granite:
+ return Til_pressure_from_entropy(density, entropy, &eos.Til_granite);
+ break;
+
+ case eos_planetary_id_Til_water:
+ return Til_pressure_from_entropy(density, entropy, &eos.Til_water);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* Hubbard & MacFarlane (1980) EoS */
+ case eos_planetary_type_HM80:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ return HM80_pressure_from_entropy(density, entropy, &eos.HM80_HHe);
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ return HM80_pressure_from_entropy(density, entropy, &eos.HM80_ice);
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ return HM80_pressure_from_entropy(density, entropy, &eos.HM80_rock);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* ANEOS EoS */
+ case eos_planetary_type_ANEOS:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_iron:
+ return ANEOS_pressure_from_entropy(density, entropy, &eos.ANEOS_iron);
+ break;
+
+ case eos_planetary_id_MANEOS_forsterite:
+ return ANEOS_pressure_from_entropy(density, entropy,
+ &eos.MANEOS_forsterite);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_pressure_from_entropy(density, entropy,
+ &eos.SESAME_iron);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ default:
+ error("Unknown material type! mat_id = %d", mat_id);
+ return 0.f;
+ }
+}
+
+/**
+ * @brief Returns the entropy given density and pressure.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param pressure The pressure \f$P\f$.
+ * @return The entropy \f$A\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
+ float density, float P, enum eos_planetary_material_id mat_id) {
+
+ const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+
+ /* Select the material base type */
+ switch (type) {
+
+ /* Tillotson EoS */
+ case eos_planetary_type_Til:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ return Til_entropy_from_pressure(density, P, &eos.Til_iron);
+ break;
+
+ case eos_planetary_id_Til_granite:
+ return Til_entropy_from_pressure(density, P, &eos.Til_granite);
+ break;
+
+ case eos_planetary_id_Til_water:
+ return Til_entropy_from_pressure(density, P, &eos.Til_water);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* Hubbard & MacFarlane (1980) EoS */
+ case eos_planetary_type_HM80:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ return HM80_entropy_from_pressure(density, P, &eos.HM80_HHe);
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ return HM80_entropy_from_pressure(density, P, &eos.HM80_ice);
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ return HM80_entropy_from_pressure(density, P, &eos.HM80_rock);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* ANEOS EoS */
+ case eos_planetary_type_ANEOS:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_iron:
+ return ANEOS_entropy_from_pressure(density, P, &eos.ANEOS_iron);
+ break;
+
+ case eos_planetary_id_MANEOS_forsterite:
+ return ANEOS_entropy_from_pressure(density, P,
+ &eos.MANEOS_forsterite);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_entropy_from_pressure(density, P, &eos.SESAME_iron);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ default:
+ error("Unknown material type! mat_id = %d", mat_id);
+ return 0.f;
+ }
+}
+
+/**
+ * @brief Returns the sound speed given density and entropy
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
+ float density, float entropy, enum eos_planetary_material_id mat_id) {
+
+ const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+
+ /* Select the material base type */
+ switch (type) {
+
+ /* Tillotson EoS */
+ case eos_planetary_type_Til:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ return Til_soundspeed_from_entropy(density, entropy, &eos.Til_iron);
+ break;
+
+ case eos_planetary_id_Til_granite:
+ return Til_soundspeed_from_entropy(density, entropy,
+ &eos.Til_granite);
+ break;
+
+ case eos_planetary_id_Til_water:
+ return Til_soundspeed_from_entropy(density, entropy, &eos.Til_water);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* Hubbard & MacFarlane (1980) EoS */
+ case eos_planetary_type_HM80:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ return HM80_soundspeed_from_entropy(density, entropy, &eos.HM80_HHe);
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ return HM80_soundspeed_from_entropy(density, entropy, &eos.HM80_ice);
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ return HM80_soundspeed_from_entropy(density, entropy, &eos.HM80_rock);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* ANEOS EoS */
+ case eos_planetary_type_ANEOS:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_iron:
+ return ANEOS_soundspeed_from_entropy(density, entropy,
+ &eos.ANEOS_iron);
+ break;
+
+ case eos_planetary_id_MANEOS_forsterite:
+ return ANEOS_soundspeed_from_entropy(density, entropy,
+ &eos.MANEOS_forsterite);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_soundspeed_from_entropy(density, entropy,
+ &eos.SESAME_iron);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ default:
+ error("Unknown material type! mat_id = %d", mat_id);
+ return 0.f;
+ }
+}
+
+/**
+ * @brief Returns the entropy given density and internal energy
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_entropy_from_internal_energy(float density, float u,
+ enum eos_planetary_material_id mat_id) {
+
+ const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+
+ /* Select the material base type */
+ switch (type) {
+
+ /* Tillotson EoS */
+ case eos_planetary_type_Til:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ return Til_entropy_from_internal_energy(density, u, &eos.Til_iron);
+ break;
+
+ case eos_planetary_id_Til_granite:
+ return Til_entropy_from_internal_energy(density, u, &eos.Til_granite);
+ break;
+
+ case eos_planetary_id_Til_water:
+ return Til_entropy_from_internal_energy(density, u, &eos.Til_water);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* Hubbard & MacFarlane (1980) EoS */
+ case eos_planetary_type_HM80:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ return HM80_entropy_from_internal_energy(density, u, &eos.HM80_HHe);
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ return HM80_entropy_from_internal_energy(density, u, &eos.HM80_ice);
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ return HM80_entropy_from_internal_energy(density, u, &eos.HM80_rock);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* ANEOS EoS */
+ case eos_planetary_type_ANEOS:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_iron:
+ return ANEOS_entropy_from_internal_energy(density, u,
+ &eos.ANEOS_iron);
+ break;
+
+ case eos_planetary_id_MANEOS_forsterite:
+ return ANEOS_entropy_from_internal_energy(density, u,
+ &eos.MANEOS_forsterite);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_entropy_from_internal_energy(density, u,
+ &eos.SESAME_iron);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ default:
+ error("Unknown material type! mat_id = %d", mat_id);
+ return 0.f;
+ }
+}
+
+/**
+ * @brief Returns the pressure given density and internal energy
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_pressure_from_internal_energy(float density, float u,
+ enum eos_planetary_material_id mat_id) {
+
+ const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+
+ /* Select the material base type */
+ switch (type) {
+
+ /* Tillotson EoS */
+ case eos_planetary_type_Til:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ return Til_pressure_from_internal_energy(density, u, &eos.Til_iron);
+ break;
+
+ case eos_planetary_id_Til_granite:
+ return Til_pressure_from_internal_energy(density, u,
+ &eos.Til_granite);
+ break;
+
+ case eos_planetary_id_Til_water:
+ return Til_pressure_from_internal_energy(density, u, &eos.Til_water);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* Hubbard & MacFarlane (1980) EoS */
+ case eos_planetary_type_HM80:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ return HM80_pressure_from_internal_energy(density, u, &eos.HM80_HHe);
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ return HM80_pressure_from_internal_energy(density, u, &eos.HM80_ice);
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ return HM80_pressure_from_internal_energy(density, u, &eos.HM80_rock);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* ANEOS EoS */
+ case eos_planetary_type_ANEOS:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_iron:
+ return ANEOS_pressure_from_internal_energy(density, u,
+ &eos.ANEOS_iron);
+ break;
+
+ case eos_planetary_id_MANEOS_forsterite:
+ return ANEOS_pressure_from_internal_energy(density, u,
+ &eos.MANEOS_forsterite);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_pressure_from_internal_energy(density, u,
+ &eos.SESAME_iron);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ default:
+ error("Unknown material type! mat_id = %d", mat_id);
+ return 0.f;
+ }
+}
+
+/**
+ * @brief Returns the internal energy given density and pressure.
+ *
+ * NOT IMPLEMENTED!
+ *
+ * @param density The density \f$\rho\f$.
+ * @param pressure The pressure \f$P\f$.
+ * @return The internal energy \f$u\f$.
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_pressure(float density, float P,
+ enum eos_planetary_material_id mat_id) {
+
+ const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+
+ /* Select the material base type */
+ switch (type) {
+
+ /* Tillotson EoS */
+ case eos_planetary_type_Til:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ return Til_internal_energy_from_pressure(density, P, &eos.Til_iron);
+ break;
+
+ case eos_planetary_id_Til_granite:
+ return Til_internal_energy_from_pressure(density, P,
+ &eos.Til_granite);
+ break;
+
+ case eos_planetary_id_Til_water:
+ return Til_internal_energy_from_pressure(density, P, &eos.Til_water);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* Hubbard & MacFarlane (1980) EoS */
+ case eos_planetary_type_HM80:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ return HM80_internal_energy_from_pressure(density, P, &eos.HM80_HHe);
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ return HM80_internal_energy_from_pressure(density, P, &eos.HM80_ice);
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ return HM80_internal_energy_from_pressure(density, P, &eos.HM80_rock);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* ANEOS EoS */
+ case eos_planetary_type_ANEOS:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_iron:
+ return ANEOS_internal_energy_from_pressure(density, P,
+ &eos.ANEOS_iron);
+ break;
+
+ case eos_planetary_id_MANEOS_forsterite:
+ return ANEOS_internal_energy_from_pressure(density, P,
+ &eos.MANEOS_forsterite);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_internal_energy_from_pressure(density, P,
+ &eos.SESAME_iron);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ default:
+ error("Unknown material type! mat_id = %d", mat_id);
+ return 0.f;
+ }
+}
+
+/**
+ * @brief Returns the sound speed given density and internal energy
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_soundspeed_from_internal_energy(float density, float u,
+ enum eos_planetary_material_id mat_id) {
+
+ const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+
+ /* Select the material base type */
+ switch (type) {
+
+ /* Tillotson EoS */
+ case eos_planetary_type_Til:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ return Til_soundspeed_from_internal_energy(density, u, &eos.Til_iron);
+ break;
+
+ case eos_planetary_id_Til_granite:
+ return Til_soundspeed_from_internal_energy(density, u,
+ &eos.Til_granite);
+ break;
+
+ case eos_planetary_id_Til_water:
+ return Til_soundspeed_from_internal_energy(density, u,
+ &eos.Til_water);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* Hubbard & MacFarlane (1980) EoS */
+ case eos_planetary_type_HM80:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ return HM80_soundspeed_from_internal_energy(density, u,
+ &eos.HM80_HHe);
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ return HM80_soundspeed_from_internal_energy(density, u,
+ &eos.HM80_ice);
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ return HM80_soundspeed_from_internal_energy(density, u,
+ &eos.HM80_rock);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* ANEOS EoS */
+ case eos_planetary_type_ANEOS:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_iron:
+ return ANEOS_soundspeed_from_internal_energy(density, u,
+ &eos.ANEOS_iron);
+ break;
+
+ case eos_planetary_id_MANEOS_forsterite:
+ return ANEOS_soundspeed_from_internal_energy(density, u,
+ &eos.MANEOS_forsterite);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_soundspeed_from_internal_energy(density, u,
+ &eos.SESAME_iron);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ default:
+ error("Unknown material type! mat_id = %d", mat_id);
+ return 0.f;
+ }
+}
+
+/**
+ * @brief Returns the sound speed given density and pressure
+ *
+ * @param density The density \f$\rho\f$
+ * @param P The pressure \f$P\f$
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
+ float density, float P, enum eos_planetary_material_id mat_id) {
+
+ const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+
+ /* Select the material base type */
+ switch (type) {
+
+ /* Tillotson EoS */
+ case eos_planetary_type_Til:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ return Til_soundspeed_from_pressure(density, P, &eos.Til_iron);
+ break;
+
+ case eos_planetary_id_Til_granite:
+ return Til_soundspeed_from_pressure(density, P, &eos.Til_granite);
+ break;
+
+ case eos_planetary_id_Til_water:
+ return Til_soundspeed_from_pressure(density, P, &eos.Til_water);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* Hubbard & MacFarlane (1980) EoS */
+ case eos_planetary_type_HM80:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ return HM80_soundspeed_from_pressure(density, P, &eos.HM80_HHe);
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ return HM80_soundspeed_from_pressure(density, P, &eos.HM80_ice);
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ return HM80_soundspeed_from_pressure(density, P, &eos.HM80_rock);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* ANEOS EoS */
+ case eos_planetary_type_ANEOS:
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_iron:
+ return ANEOS_soundspeed_from_pressure(density, P, &eos.ANEOS_iron);
+ break;
+
+ case eos_planetary_id_MANEOS_forsterite:
+ return ANEOS_soundspeed_from_pressure(density, P,
+ &eos.MANEOS_forsterite);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
+
+ /* Select the material */
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_soundspeed_from_pressure(density, P, &eos.SESAME_iron);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ return 0.f;
+ };
+ break;
+
+ default:
+ error("Unknown material type! mat_id = %d", mat_id);
+ return 0.f;
+ }
+}
+
+/**
+ * @brief Initialize the eos parameters
+ *
+ * @param e The #eos_parameters
+ * @param params The parsed parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_init(
+ struct eos_parameters *e, const struct phys_const *phys_const,
+ const struct unit_system *us, const struct swift_params *params) {
+
+ // Set the parameters and material IDs, load tables, etc. for each material
+ // Tillotson
+ set_Til_iron(&e->Til_iron, eos_planetary_id_Til_iron);
+ set_Til_granite(&e->Til_granite, eos_planetary_id_Til_granite);
+ set_Til_water(&e->Til_water, eos_planetary_id_Til_water);
+
+ // Hubbard & MacFarlane (1980)
+ set_HM80_HHe(&e->HM80_HHe, eos_planetary_id_HM80_HHe);
+ set_HM80_ice(&e->HM80_ice, eos_planetary_id_HM80_ice);
+ set_HM80_rock(&e->HM80_rock, eos_planetary_id_HM80_rock);
+
+ load_HM80_table(&e->HM80_HHe, HM80_HHe_table_file);
+ load_HM80_table(&e->HM80_ice, HM80_ice_table_file);
+ load_HM80_table(&e->HM80_rock, HM80_rock_table_file);
+
+ // ANEOS
+ set_ANEOS_iron(&e->ANEOS_iron, eos_planetary_id_ANEOS_iron);
+ set_MANEOS_forsterite(&e->MANEOS_forsterite,
+ eos_planetary_id_MANEOS_forsterite);
+
+ // SESAME
+ set_SESAME_iron(&e->SESAME_iron, eos_planetary_id_SESAME_iron);
+
+ // Convert to internal units
+ // Tillotson
+ convert_units_Til(&e->Til_iron, us);
+ convert_units_Til(&e->Til_granite, us);
+ convert_units_Til(&e->Til_water, us);
+
+ // Hubbard & MacFarlane (1980)
+ convert_units_HM80(&e->HM80_HHe, us);
+ convert_units_HM80(&e->HM80_ice, us);
+ convert_units_HM80(&e->HM80_rock, us);
+
+ // ANEOS
+ convert_units_ANEOS(&e->ANEOS_iron, us);
+ convert_units_ANEOS(&e->MANEOS_forsterite, us);
+
+ // SESAME
+ convert_units_SESAME(&e->SESAME_iron, us);
+}
+
+/**
+ * @brief Print the equation of state
+ *
+ * @param e The #eos_parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_print(
+ const struct eos_parameters *e) {
+
+ message("Equation of state: Planetary.");
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Write equation of state information to the snapshot
+ *
+ * @param h_grpsph The HDF5 group in which to write
+ * @param e The #eos_parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_print_snapshot(
+ hid_t h_grpsph, const struct eos_parameters *e) {
+
+ io_write_attribute_s(h_grpsph, "Equation of state", "Planetary");
+}
+#endif
+
+#endif /* SWIFT_PLANETARY_EQUATION_OF_STATE_H */
diff --git a/src/equation_of_state/planetary/hm80.h b/src/equation_of_state/planetary/hm80.h
new file mode 100644
index 0000000000000000000000000000000000000000..b0669c129b2b4d642ae27e58f2f2bec221e3df9b
--- /dev/null
+++ b/src/equation_of_state/planetary/hm80.h
@@ -0,0 +1,306 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * 2018 Jacob Kegerreis (jacob.kegerreis@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_HUBBARD_MACFARLANE_EQUATION_OF_STATE_H
+#define SWIFT_HUBBARD_MACFARLANE_EQUATION_OF_STATE_H
+
+/**
+ * @file equation_of_state/planetary/hm80.h
+ *
+ * Contains the Hubbard & MacFarlane (1980) Uranus/Neptune EOS functions for
+ * equation_of_state/planetary/equation_of_state.h
+ *
+ */
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "common_io.h"
+#include "equation_of_state.h"
+#include "inline.h"
+#include "physical_constants.h"
+#include "units.h"
+
+// Hubbard & MacFarlane (1980) parameters
+struct HM80_params {
+ float **table_P_rho_u;
+ int num_rho, num_u;
+ float log_rho_min, log_rho_max, log_rho_step, inv_log_rho_step, log_u_min,
+ log_u_max, log_u_step, inv_log_u_step, bulk_mod;
+ enum eos_planetary_material_id mat_id;
+};
+
+// Table file names
+/// to be read in from the parameter file instead once finished testing...
+#define HM80_HHe_table_file "/gpfs/data/dc-kege1/gihr_data/P_rho_u_HHe.txt"
+#define HM80_ice_table_file "/gpfs/data/dc-kege1/gihr_data/P_rho_u_ice.txt"
+#define HM80_rock_table_file "/gpfs/data/dc-kege1/gihr_data/P_rho_u_roc.txt"
+
+// Parameter values for each material (cgs units)
+INLINE static void set_HM80_HHe(struct HM80_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ mat->mat_id = mat_id;
+ mat->num_rho = 100;
+ mat->num_u = 100;
+ mat->log_rho_min = -9.2103404f;
+ mat->log_rho_max = 1.6094379f;
+ mat->log_rho_step = 0.1092907f;
+ mat->log_u_min = 9.2103404f;
+ mat->log_u_max = 22.3327037f;
+ mat->log_u_step = 0.1325491f;
+ mat->bulk_mod = 0;
+
+ mat->inv_log_rho_step = 1.f / mat->log_rho_step;
+ mat->inv_log_u_step = 1.f / mat->log_u_step;
+}
+INLINE static void set_HM80_ice(struct HM80_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ mat->mat_id = mat_id;
+ mat->num_rho = 200;
+ mat->num_u = 200;
+ mat->log_rho_min = -6.9077553f;
+ mat->log_rho_max = 2.7080502f;
+ mat->log_rho_step = 0.0483206f;
+ mat->log_u_min = 6.9077553f;
+ mat->log_u_max = 22.3327037f;
+ mat->log_u_step = 0.0775123f;
+ mat->bulk_mod = 2.0e10f;
+
+ mat->inv_log_rho_step = 1.f / mat->log_rho_step;
+ mat->inv_log_u_step = 1.f / mat->log_u_step;
+}
+INLINE static void set_HM80_rock(struct HM80_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ mat->mat_id = mat_id;
+ mat->num_rho = 100;
+ mat->num_u = 100;
+ mat->log_rho_min = -6.9077553f;
+ mat->log_rho_max = 2.9957323f;
+ mat->log_rho_step = 0.1000352f;
+ mat->log_u_min = 9.2103404f;
+ mat->log_u_max = 20.7232658f;
+ mat->log_u_step = 0.1162922f;
+ mat->bulk_mod = 3.49e11f;
+
+ mat->inv_log_rho_step = 1.f / mat->log_rho_step;
+ mat->inv_log_u_step = 1.f / mat->log_u_step;
+}
+
+// Read the table from file
+INLINE static void load_HM80_table(struct HM80_params *mat, char *table_file) {
+ // Allocate table memory
+ mat->table_P_rho_u = (float **)malloc(mat->num_rho * sizeof(float *));
+ for (int i = 0; i < mat->num_rho; i++) {
+ mat->table_P_rho_u[i] = (float *)malloc(mat->num_u * sizeof(float));
+ }
+
+ // Load table contents from file
+ FILE *f = fopen(table_file, "r");
+ int c;
+ for (int i = 0; i < mat->num_rho; i++) {
+ for (int j = 0; j < mat->num_u; j++) {
+ c = fscanf(f, "%f", &mat->table_P_rho_u[i][j]);
+ if (c != 1) {
+ error("Failed to read EOS table");
+ }
+ }
+ }
+ fclose(f);
+}
+
+// Convert from cgs to internal units
+INLINE static void convert_units_HM80(struct HM80_params *mat,
+ const struct unit_system *us) {
+ const float Mbar_to_Ba = 1e12f; // Convert Megabar to Barye
+ const float J_kg_to_erg_g = 1e4f; // Convert J/kg to erg/g
+
+ // Table densities in cgs
+ mat->log_rho_min -= logf(units_cgs_conversion_factor(us, UNIT_CONV_DENSITY));
+ mat->log_rho_max -= logf(units_cgs_conversion_factor(us, UNIT_CONV_DENSITY));
+
+ // Table energies in SI
+ mat->log_u_min +=
+ logf(J_kg_to_erg_g /
+ units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
+ mat->log_u_max +=
+ logf(J_kg_to_erg_g /
+ units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
+
+ // Table Pressures in Mbar
+ for (int i = 0; i < mat->num_rho; i++) {
+ for (int j = 0; j < mat->num_u; j++) {
+ mat->table_P_rho_u[i][j] *=
+ Mbar_to_Ba / units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
+ }
+ }
+
+ mat->bulk_mod /= units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
+}
+
+// gas_internal_energy_from_entropy
+INLINE static float HM80_internal_energy_from_entropy(
+ float density, float entropy, const struct HM80_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_pressure_from_entropy
+INLINE static float HM80_pressure_from_entropy(float density, float entropy,
+ const struct HM80_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_entropy_from_pressure
+INLINE static float HM80_entropy_from_pressure(float density, float pressure,
+ const struct HM80_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_soundspeed_from_entropy
+INLINE static float HM80_soundspeed_from_entropy(
+ float density, float entropy, const struct HM80_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_entropy_from_internal_energy
+INLINE static float HM80_entropy_from_internal_energy(
+ float density, float u, const struct HM80_params *mat) {
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_pressure_from_internal_energy
+INLINE static float HM80_pressure_from_internal_energy(
+ float density, float u, const struct HM80_params *mat) {
+
+ float P;
+
+ if (u <= 0) {
+ return 0.f;
+ }
+
+ int rho_idx, u_idx;
+ float intp_rho, intp_u;
+ const float log_rho = logf(density);
+ const float log_u = logf(u);
+
+ // 2D interpolation (linear in log(rho), log(u)) to find P(rho, u)
+ rho_idx = floorf((log_rho - mat->log_rho_min) * mat->inv_log_rho_step);
+ u_idx = floorf((log_u - mat->log_u_min) * mat->inv_log_u_step);
+
+ intp_rho = (log_rho - mat->log_rho_min - rho_idx * mat->log_rho_step) *
+ mat->inv_log_rho_step;
+ intp_u =
+ (log_u - mat->log_u_min - u_idx * mat->log_u_step) * mat->inv_log_u_step;
+
+ // Return zero pressure if below the table minimum/a
+ // Extrapolate the pressure for low densities
+ if (rho_idx < 0) { // Too-low rho
+ P = expf(logf((1 - intp_u) * mat->table_P_rho_u[0][u_idx] +
+ intp_u * mat->table_P_rho_u[0][u_idx + 1]) +
+ log_rho - mat->log_rho_min);
+ if (u_idx < 0) { // and too-low u
+ P = 0;
+ }
+ } else if (u_idx < 0) { // Too-low u
+ P = 0;
+ }
+ // Return an edge value if above the table maximum/a
+ else if (rho_idx >= mat->num_rho - 1) { // Too-high rho
+ if (u_idx >= mat->num_u - 1) { // and too-high u
+ P = mat->table_P_rho_u[mat->num_rho - 1][mat->num_u - 1];
+ } else {
+ P = mat->table_P_rho_u[mat->num_rho - 1][u_idx];
+ }
+ } else if (u_idx >= mat->num_u - 1) { // Too-high u
+ P = mat->table_P_rho_u[rho_idx][mat->num_u - 1];
+ }
+ // Normal interpolation within the table
+ else {
+ P = (1.f - intp_rho) *
+ ((1.f - intp_u) * mat->table_P_rho_u[rho_idx][u_idx] +
+ intp_u * mat->table_P_rho_u[rho_idx][u_idx + 1]) +
+ intp_rho * ((1 - intp_u) * mat->table_P_rho_u[rho_idx + 1][u_idx] +
+ intp_u * mat->table_P_rho_u[rho_idx + 1][u_idx + 1]);
+ }
+
+ return P;
+}
+
+// gas_internal_energy_from_pressure
+INLINE static float HM80_internal_energy_from_pressure(
+ float density, float P, const struct HM80_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_soundspeed_from_internal_energy
+INLINE static float HM80_soundspeed_from_internal_energy(
+ float density, float u, const struct HM80_params *mat) {
+
+ float c, P;
+
+ // Bulk modulus
+ if (mat->bulk_mod != 0) {
+ c = sqrtf(mat->bulk_mod / density);
+ }
+ // Ideal gas
+ else {
+ P = HM80_pressure_from_internal_energy(density, u, mat);
+ c = sqrtf(hydro_gamma * P / density);
+ }
+
+ return c;
+}
+
+// gas_soundspeed_from_pressure
+INLINE static float HM80_soundspeed_from_pressure(
+ float density, float P, const struct HM80_params *mat) {
+
+ float c;
+
+ // Bulk modulus
+ if (mat->bulk_mod != 0) {
+ c = sqrtf(mat->bulk_mod / density);
+ }
+ // Ideal gas
+ else {
+ c = sqrtf(hydro_gamma * P / density);
+ }
+
+ return c;
+}
+
+#endif /* SWIFT_HUBBARD_MACFARLANE_EQUATION_OF_STATE_H */
diff --git a/src/equation_of_state/planetary/sesame.h b/src/equation_of_state/planetary/sesame.h
new file mode 100644
index 0000000000000000000000000000000000000000..76574c2ad00282a82649705cd8a2b5a1b428d867
--- /dev/null
+++ b/src/equation_of_state/planetary/sesame.h
@@ -0,0 +1,140 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * 2018 Jacob Kegerreis (jacob.kegerreis@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_SESAME_EQUATION_OF_STATE_H
+#define SWIFT_SESAME_EQUATION_OF_STATE_H
+
+/**
+ * @file equation_of_state/planetary/sesame.h
+ *
+ * Contains the SESAME EOS functions for
+ * equation_of_state/planetary/equation_of_state.h
+ *
+ * WORK IN PROGRESS!
+ *
+ */
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "common_io.h"
+#include "equation_of_state.h"
+#include "inline.h"
+#include "physical_constants.h"
+#include "units.h"
+
+// SESAME parameters
+struct SESAME_params {
+ enum eos_planetary_material_id mat_id;
+};
+
+// Parameter values for each material (cgs units)
+INLINE static void set_SESAME_iron(struct SESAME_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ mat->mat_id = mat_id;
+}
+
+// Convert from cgs to internal units
+INLINE static void convert_units_SESAME(struct SESAME_params *mat,
+ const struct unit_system *us) {}
+
+// gas_internal_energy_from_entropy
+INLINE static float SESAME_internal_energy_from_entropy(
+ float density, float entropy, const struct SESAME_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_pressure_from_entropy
+INLINE static float SESAME_pressure_from_entropy(
+ float density, float entropy, const struct SESAME_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_entropy_from_pressure
+INLINE static float SESAME_entropy_from_pressure(
+ float density, float pressure, const struct SESAME_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_soundspeed_from_entropy
+INLINE static float SESAME_soundspeed_from_entropy(
+ float density, float entropy, const struct SESAME_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_entropy_from_internal_energy
+INLINE static float SESAME_entropy_from_internal_energy(
+ float density, float u, const struct SESAME_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_pressure_from_internal_energy
+INLINE static float SESAME_pressure_from_internal_energy(
+ float density, float u, const struct SESAME_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_internal_energy_from_pressure
+INLINE static float SESAME_internal_energy_from_pressure(
+ float density, float P, const struct SESAME_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_soundspeed_from_internal_energy
+INLINE static float SESAME_soundspeed_from_internal_energy(
+ float density, float u, const struct SESAME_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_soundspeed_from_pressure
+INLINE static float SESAME_soundspeed_from_pressure(
+ float density, float P, const struct SESAME_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+#endif /* SWIFT_SESAME_EQUATION_OF_STATE_H */
diff --git a/src/equation_of_state/planetary/tillotson.h b/src/equation_of_state/planetary/tillotson.h
new file mode 100644
index 0000000000000000000000000000000000000000..ed280079b040776dfc20480c50fbc9a3b6ff4470
--- /dev/null
+++ b/src/equation_of_state/planetary/tillotson.h
@@ -0,0 +1,282 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * 2018 Jacob Kegerreis (jacob.kegerreis@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_TILLOTSON_EQUATION_OF_STATE_H
+#define SWIFT_TILLOTSON_EQUATION_OF_STATE_H
+
+/**
+ * @file equation_of_state/planetary/tillotson.h
+ *
+ * Contains the Tillotson EOS functions for
+ * equation_of_state/planetary/equation_of_state.h
+ *
+ */
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "common_io.h"
+#include "equation_of_state.h"
+#include "inline.h"
+#include "physical_constants.h"
+#include "units.h"
+
+// Tillotson parameters
+struct Til_params {
+ float rho_0, a, b, A, B, E_0, E_iv, E_cv, alpha, beta, eta_min, P_min;
+ enum eos_planetary_material_id mat_id;
+};
+
+// Parameter values for each material (cgs units)
+INLINE static void set_Til_iron(struct Til_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ mat->mat_id = mat_id;
+ mat->rho_0 = 7.800f;
+ mat->a = 0.5f;
+ mat->b = 1.5f;
+ mat->A = 1.28e12f;
+ mat->B = 1.05e12f;
+ mat->E_0 = 9.5e10f;
+ mat->E_iv = 2.4e10f;
+ mat->E_cv = 8.67e10f;
+ mat->alpha = 5.0f;
+ mat->beta = 5.0f;
+ mat->eta_min = 0.0f;
+ mat->P_min = 0.0f;
+}
+INLINE static void set_Til_granite(struct Til_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ mat->mat_id = mat_id;
+ mat->rho_0 = 2.680f;
+ mat->a = 0.5f;
+ mat->b = 1.3f;
+ mat->A = 1.8e11f;
+ mat->B = 1.8e11f;
+ mat->E_0 = 1.6e11f;
+ mat->E_iv = 3.5e10f;
+ mat->E_cv = 1.8e11f;
+ mat->alpha = 5.0f;
+ mat->beta = 5.0f;
+ mat->eta_min = 0.0f;
+ mat->P_min = 0.0f;
+}
+INLINE static void set_Til_water(struct Til_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ mat->mat_id = mat_id;
+ mat->rho_0 = 0.998f;
+ mat->a = 0.7f;
+ mat->b = 0.15f;
+ mat->A = 2.18e10f;
+ mat->B = 1.325e11f;
+ mat->E_0 = 7.0e10f;
+ mat->E_iv = 4.19e9f;
+ mat->E_cv = 2.69e10f;
+ mat->alpha = 10.0f;
+ mat->beta = 5.0f;
+ mat->eta_min = 0.915f;
+ mat->P_min = 0.0f;
+}
+
+// Convert from cgs to internal units
+INLINE static void convert_units_Til(struct Til_params *mat,
+ const struct unit_system *us) {
+
+ mat->rho_0 /= units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+ mat->A /= units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
+ mat->B /= units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
+ mat->E_0 /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ mat->E_iv /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ mat->E_cv /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ mat->P_min /= units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
+}
+
+// gas_internal_energy_from_entropy
+INLINE static float Til_internal_energy_from_entropy(
+ float density, float entropy, const struct Til_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_pressure_from_entropy
+INLINE static float Til_pressure_from_entropy(float density, float entropy,
+ const struct Til_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_entropy_from_pressure
+INLINE static float Til_entropy_from_pressure(float density, float pressure,
+ const struct Til_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_soundspeed_from_entropy
+INLINE static float Til_soundspeed_from_entropy(float density, float entropy,
+ const struct Til_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_entropy_from_internal_energy
+INLINE static float Til_entropy_from_internal_energy(
+ float density, float u, const struct Til_params *mat) {
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_pressure_from_internal_energy
+INLINE static float Til_pressure_from_internal_energy(
+ float density, float u, const struct Til_params *mat) {
+
+ const float eta = density / mat->rho_0;
+ const float mu = eta - 1.f;
+ const float nu = 1.f / eta - 1.f;
+ float P_c, P_e, P;
+
+ // Condensed or cold
+ if (eta < mat->eta_min) {
+ P_c = 0.f;
+ } else {
+ P_c = (mat->a + mat->b / (u / (mat->E_0 * eta * eta) + 1.f)) * density * u +
+ mat->A * mu + mat->B * mu * mu;
+ }
+ // Expanded and hot
+ P_e = mat->a * density * u +
+ (mat->b * density * u / (u / (mat->E_0 * eta * eta) + 1.f) +
+ mat->A * mu * expf(-mat->beta * nu)) *
+ expf(-mat->alpha * nu * nu);
+
+ // Condensed or cold state
+ if ((1.f < eta) || (u < mat->E_iv)) {
+ P = P_c;
+ }
+ // Expanded and hot state
+ else if ((eta < 1.f) && (mat->E_cv < u)) {
+ P = P_e;
+ }
+ // Hybrid state
+ else {
+ P = ((u - mat->E_iv) * P_e + (mat->E_cv - u) * P_c) /
+ (mat->E_cv - mat->E_iv);
+ }
+
+ // Minimum pressure
+ if (P < mat->P_min) {
+ P = mat->P_min;
+ }
+
+ return P;
+}
+
+// gas_internal_energy_from_pressure
+INLINE static float Til_internal_energy_from_pressure(
+ float density, float P, const struct Til_params *mat) {
+
+ error("This EOS function is not yet implemented!");
+
+ return 0;
+}
+
+// gas_soundspeed_from_internal_energy
+INLINE static float Til_soundspeed_from_internal_energy(
+ float density, float u, const struct Til_params *mat) {
+
+ // const float eta = density / mat->rho_0;
+ // const float mu = eta - 1.f;
+ // const float nu = 1.f/eta - 1.f;
+ // float P_c, P_e, P, c_c, c_e, c;
+ //
+ // // Condensed or cold
+ // if (eta < mat->eta_min) {
+ // P_c = 0.f;
+ // }
+ // else {
+ // P_c = (mat->a + mat->b / (u / (mat->E_0 * eta*eta) + 1.f)) * density
+ // * u
+ // + mat->A * mu + mat->B * mu*mu;
+ // }
+ // c_c = mat->a*u + mat->b*u / ((u / (mat->E_0*eta*eta)+1.f) *
+ // (u / (mat->E_0*eta*eta)+1.f)) *
+ // (3.f*(u / (mat->E_0*eta*eta)+1.f) - 2.f) +
+ // (mat->A + 2.f*mat->B*mu) / mat->rho_0 + P_c / (rho*rho) *
+ // (mat->a*rho + mat->b*rho / ((u / (mat->E_0*eta*eta)+1.f) *
+ // (u / (mat->E_0*eta*eta)+1.f)));
+ //
+ // c_c = max(c_c, mat->A / mat->rho_0);
+ //
+ // // Expanded and hot
+ // P_e = mat->a*density*u + (
+ // mat->b * density * u / (u / (mat->E_0 * eta*eta) + 1.f)
+ // + mat->A*mu * expf(-mat->beta * nu)
+ // ) * expf(-mat->alpha * nu*nu);
+ //
+ // c_e = (mat->a + mat->b / (u / (mat->E_0*eta*eta)+1.f) *
+ // expf(-mat->beta*((1.f - eta)/eta)*((1.f - eta)/eta))
+ // + 1.f)*P_e/rho + mat->A/mat->rho_0
+ // *expf(-(mat->alpha*((1.f - eta)/eta)+mat->beta *
+ // ((1.f - eta)/eta)*((1.f - eta)/eta)))*(1.f+mu/(eta*eta)
+ // *(mat->alpha+2.f*mat->beta*((1.f - eta)/eta)-eta)) +
+ // mat->b*rho*u/((u / (mat->E_0*eta*eta)+1.f)*
+ // (u / (mat->E_0*eta*eta)+1.f)*eta*eta)*
+ // expf(-mat->beta*((1.f - eta)/eta)*((1.f - eta)/eta))*
+ // (2.f*mat->beta*((1.f - eta)/eta)*(u / (mat->E_0*eta*eta)+1.f) /
+ // mat->rho_0 + 1.f/(mat->E_0*rho)*(2.f*u-P_e/rho));
+ //
+ // // Condensed or cold state
+ // if ((1.f < eta) || (u < mat->E_iv)) {
+ // c = c_c;
+ // }
+ // // Expanded and hot state
+ // else if ((eta < 1.f) && (mat->E_cv < u)) {
+ // c = c_e;
+ // }
+ // // Hybrid state
+ // else {
+ // c = ((u - mat->E_iv)*c_e + (mat->E_cv - u)*c_c) /
+ // (mat->E_cv - mat->E_iv);
+ //
+ // c = max(c_c, mat->A / mat->rho0);
+ // }
+ float c = sqrtf(mat->A / mat->rho_0);
+
+ return c;
+}
+
+// gas_soundspeed_from_pressure
+INLINE static float Til_soundspeed_from_pressure(float density, float P,
+ const struct Til_params *mat) {
+
+ float c = sqrtf(mat->A / mat->rho_0);
+
+ return c;
+}
+
+#endif /* SWIFT_TILLOTSON_EQUATION_OF_STATE_H */
diff --git a/src/hydro.h b/src/hydro.h
index 78ae7d178ff93b01fa2eeebe09f34b718ce9a826..0b6438471ec799f9b44e01b48dd528de9f57e379 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -54,6 +54,10 @@
#include "./hydro/Shadowswift/hydro_iact.h"
#define SPH_IMPLEMENTATION \
"Shadowfax moving mesh (Vandenbroucke and De Rijcke 2016)"
+#elif defined(MINIMAL_MULTI_MAT_SPH)
+#include "./hydro/MinimalMultiMat/hydro.h"
+#include "./hydro/MinimalMultiMat/hydro_iact.h"
+#define SPH_IMPLEMENTATION "Minimal version of SPH with multiple materials"
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/hydro/MinimalMultiMat/hydro.h b/src/hydro/MinimalMultiMat/hydro.h
new file mode 100644
index 0000000000000000000000000000000000000000..5383ffda8fe67a591691766e4150e75a9dbd4cb0
--- /dev/null
+++ b/src/hydro/MinimalMultiMat/hydro.h
@@ -0,0 +1,634 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2018 Jacob Kegerreis (jacob.kegerreis@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_MINIMAL_MULTI_MAT_HYDRO_H
+#define SWIFT_MINIMAL_MULTI_MAT_HYDRO_H
+
+/**
+ * @file MinimalMultiMat/hydro.h
+ * @brief Minimal conservative implementation of SPH (Non-neighbour loop
+ * equations) with multiple materials.
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
+ * Physics, 2012, Volume 231, Issue 3, pp. 759-794.
+ */
+
+#include "adiabatic_index.h"
+#include "approx_math.h"
+#include "cosmology.h"
+#include "dimension.h"
+#include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "hydro_space.h"
+#include "kernel_hydro.h"
+#include "minmax.h"
+
+/**
+ * @brief Returns the comoving internal energy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy(const struct part *restrict p) {
+
+ return p->u;
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable and converts it to
+ * physical coordinates.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return p->u * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Returns the comoving pressure of a particle
+ *
+ * Computes the pressure based on the particle's properties.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
+ const struct part *restrict p) {
+
+ return gas_pressure_from_internal_energy(p->rho, p->u, p->mat_id);
+}
+
+/**
+ * @brief Returns the physical pressure of a particle
+ *
+ * Computes the pressure based on the particle's properties and
+ * convert it to physical coordinates.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ return cosmo->a_factor_pressure *
+ gas_pressure_from_internal_energy(p->rho, p->u, p->mat_id);
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
+ const struct part *restrict p) {
+
+ return gas_entropy_from_internal_energy(p->rho, p->u, p->mat_id);
+}
+
+/**
+ * @brief Returns the physical entropy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable and converts it to
+ * physical coordinates.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ /* Note: no cosmological conversion required here with our choice of
+ * coordinates. */
+ return gas_entropy_from_internal_energy(p->rho, p->u, p->mat_id);
+}
+
+/**
+ * @brief Returns the comoving sound speed of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_soundspeed(const struct part *restrict p) {
+
+ return p->force.soundspeed;
+}
+
+/**
+ * @brief Returns the physical sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_soundspeed(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return cosmo->a_factor_sound_speed * p->force.soundspeed;
+}
+
+/**
+ * @brief Returns the comoving density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
+ const struct part *restrict p) {
+
+ return p->rho;
+}
+
+/**
+ * @brief Returns the comoving density of a particle.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ return cosmo->a3_inv * p->rho;
+}
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+ const struct part *restrict p) {
+
+ return p->mass;
+}
+
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part *restrict p, float m) {
+
+ p->mass = m;
+}
+
+/**
+ * @brief Returns the velocities drifted to the current time of a particle.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle.
+ * @param dt_kick_hydro The time (for hydro accelerations) since the last kick.
+ * @param dt_kick_grav The time (for gravity accelerations) since the last kick.
+ * @param v (return) The velocities at the current time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
+ const struct part *restrict p, const struct xpart *xp, float dt_kick_hydro,
+ float dt_kick_grav, float v[3]) {
+
+ v[0] = xp->v_full[0] + p->a_hydro[0] * dt_kick_hydro +
+ xp->a_grav[0] * dt_kick_grav;
+ v[1] = xp->v_full[1] + p->a_hydro[1] * dt_kick_hydro +
+ xp->a_grav[1] * dt_kick_grav;
+ v[2] = xp->v_full[2] + p->a_hydro[2] * dt_kick_hydro +
+ xp->a_grav[2] * dt_kick_grav;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
+ const struct part *restrict p) {
+
+ return p->u_dt;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest.
+ * @param du_dt The new time derivative of the internal energy.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
+ struct part *restrict p, float du_dt) {
+
+ p->u_dt = du_dt;
+}
+/**
+ * @brief Computes the hydro time-step of a given particle
+ *
+ * This function returns the time-step of a particle given its hydro-dynamical
+ * state. A typical time-step calculation would be the use of the CFL condition.
+ *
+ * @param p Pointer to the particle data
+ * @param xp Pointer to the extended particle data
+ * @param hydro_properties The SPH parameters
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
+ const struct part *restrict p, const struct xpart *restrict xp,
+ const struct hydro_props *restrict hydro_properties,
+ const struct cosmology *restrict cosmo) {
+
+ const float CFL_condition = hydro_properties->CFL_condition;
+
+ /* CFL condition */
+ const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
+ (cosmo->a_factor_sound_speed * p->force.v_sig);
+
+ return dt_cfl;
+}
+
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part *p, float dt) {}
+
+/**
+ * @brief Prepares a particle for the density calculation.
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the various density loop over neighbours. Typically, all fields of the
+ * density sub-structure of a particle get zeroed in here.
+ *
+ * @param p The particle to act upon
+ * @param hs #hydro_space containing hydro specific space information.
+ */
+__attribute__((always_inline)) INLINE static void hydro_init_part(
+ struct part *restrict p, const struct hydro_space *hs) {
+
+ p->density.wcount = 0.f;
+ p->density.wcount_dh = 0.f;
+ p->rho = 0.f;
+ p->density.rho_dh = 0.f;
+}
+
+/**
+ * @brief Finishes the density calculation.
+ *
+ * Multiplies the density and number of neighbours by the appropiate constants
+ * and add the self-contribution term.
+ * Additional quantities such as velocity gradients will also get the final
+ * terms added to them here.
+ *
+ * Also adds/multiplies the cosmological terms if need be.
+ *
+ * @param p The particle to act upon
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+ struct part *restrict p, const struct cosmology *cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+ /* Final operation on the density (add self-contribution). */
+ p->rho += p->mass * kernel_root;
+ p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
+ p->density.wcount += kernel_root;
+ p->density.wcount_dh -= hydro_dimension * kernel_root;
+
+ /* Finish the calculation by inserting the missing h-factors */
+ p->rho *= h_inv_dim;
+ p->density.rho_dh *= h_inv_dim_plus_one;
+ p->density.wcount *= h_inv_dim;
+ p->density.wcount_dh *= h_inv_dim_plus_one;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * In the desperate case where a particle has no neighbours (likely because
+ * of the h_max ceiling), set the particle fields to something sensible to avoid
+ * NaNs in the next calculations.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+ /* Re-set problematic values */
+ p->rho = p->mass * kernel_root * h_inv_dim;
+ p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+ p->density.rho_dh = 0.f;
+ p->density.wcount_dh = 0.f;
+}
+
+/**
+ * @brief Prepare a particle for the force calculation.
+ *
+ * This function is called in the ghost task to convert some quantities coming
+ * from the density loop over neighbours into quantities ready to be used in the
+ * force loop over neighbours. Quantities are typically read from the density
+ * sub-structure and written to the force sub-structure.
+ * Examples of calculations done here include the calculation of viscosity term
+ * constants, thermal conduction terms, hydro conversions, etc.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ /* Compute the pressure */
+ const float pressure =
+ gas_pressure_from_internal_energy(p->rho, p->u, p->mat_id);
+
+ /* Compute the sound speed */
+ const float soundspeed =
+ gas_soundspeed_from_pressure(p->rho, pressure, p->mat_id);
+
+ /* Compute the "grad h" term */
+ const float rho_inv = 1.f / p->rho;
+ const float grad_h_term =
+ 1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
+
+ /* Update variables. */
+ p->force.f = grad_h_term;
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * Resets all hydro acceleration and time derivative fields in preparation
+ * for the sums taking place in the various force tasks.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+ struct part *restrict p) {
+
+ /* Reset the acceleration. */
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
+ /* Reset the time derivatives. */
+ p->u_dt = 0.0f;
+ p->force.h_dt = 0.0f;
+ p->force.v_sig = 0.0f;
+}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param p The particle.
+ * @param xp The extended data of this particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
+ struct part *restrict p, const struct xpart *restrict xp) {
+
+ /* Re-set the predicted velocities */
+ p->v[0] = xp->v_full[0];
+ p->v[1] = xp->v_full[1];
+ p->v[2] = xp->v_full[2];
+
+ /* Re-set the entropy */
+ p->u = xp->u_full;
+}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * Additional hydrodynamic quantites are drifted forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * Note the different time-step sizes used for the different quantities as they
+ * include cosmological factors.
+ *
+ * @param p The particle.
+ * @param xp The extended data of the particle.
+ * @param dt_drift The drift time-step for positions.
+ * @param dt_therm The drift time-step for thermal quantities.
+ */
+__attribute__((always_inline)) INLINE static void hydro_predict_extra(
+ struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
+ float dt_therm) {
+
+ const float h_inv = 1.f / p->h;
+
+ /* Predict smoothing length */
+ const float w1 = p->force.h_dt * h_inv * dt_drift;
+ if (fabsf(w1) < 0.2f)
+ p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+ else
+ p->h *= expf(w1);
+
+ /* Predict density */
+ const float w2 = -hydro_dimension * w1;
+ if (fabsf(w2) < 0.2f)
+ p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
+ else
+ p->rho *= expf(w2);
+
+ /* Predict the internal energy */
+ p->u += p->u_dt * dt_therm;
+
+ /* Compute the new pressure */
+ const float pressure =
+ gas_pressure_from_internal_energy(p->rho, p->u, p->mat_id);
+
+ /* Compute the new sound speed */
+ const float soundspeed =
+ gas_soundspeed_from_pressure(p->rho, pressure, p->mat_id);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+}
+
+/**
+ * @brief Finishes the force calculation.
+ *
+ * Multiplies the force and accelerations by the appropiate constants
+ * and add the self-contribution term. In most cases, there is little
+ * to do here.
+ *
+ * Cosmological terms are also added/multiplied here.
+ *
+ * @param p The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+ struct part *restrict p, const struct cosmology *cosmo) {
+
+ p->force.h_dt *= p->h * hydro_dimension_inv;
+}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * Additional hydrodynamic quantites are kicked forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * @param p The particle to act upon.
+ * @param xp The particle extended data to act upon.
+ * @param dt_therm The time-step for this kick (for thermodynamic quantities).
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+ struct part *restrict p, struct xpart *restrict xp, float dt_therm,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
+
+ /* Do not decrease the energy by more than a factor of 2*/
+ if (dt_therm > 0. && p->u_dt * dt_therm < -0.5f * xp->u_full) {
+ p->u_dt = -0.5f * xp->u_full / dt_therm;
+ }
+ xp->u_full += p->u_dt * dt_therm;
+
+ /* Apply the minimal energy limit */
+ const float min_energy =
+ hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
+ if (xp->u_full < min_energy) {
+ xp->u_full = min_energy;
+ p->u_dt = 0.f;
+ }
+
+ /* Compute the pressure */
+ const float pressure =
+ gas_pressure_from_internal_energy(p->rho, xp->u_full, p->mat_id);
+
+ /* Compute the sound speed */
+ const float soundspeed =
+ gas_soundspeed_from_internal_energy(p->rho, p->u, p->mat_id);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+}
+
+/**
+ * @brief Converts hydro quantity of a particle at the start of a run
+ *
+ * This function is called once at the end of the engine_init_particle()
+ * routine (at the start of a calculation) after the densities of
+ * particles have been computed.
+ * This can be used to convert internal energy into entropy for instance.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle to act upon
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ /* Compute the pressure */
+ const float pressure =
+ gas_pressure_from_internal_energy(p->rho, p->u, p->mat_id);
+
+ /* Compute the sound speed */
+ const float soundspeed =
+ gas_soundspeed_from_internal_energy(p->rho, p->u, p->mat_id);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+}
+
+/**
+ * @brief Initialises the particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions or assignments between the particle
+ * and extended particle fields.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+ struct part *restrict p, struct xpart *restrict xp) {
+
+ p->time_bin = 0;
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
+ xp->a_grav[0] = 0.f;
+ xp->a_grav[1] = 0.f;
+ xp->a_grav[2] = 0.f;
+ xp->u_full = p->u;
+
+ hydro_reset_acceleration(p);
+ hydro_init_part(p, NULL);
+}
+
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+ p->u = u_init;
+}
+
+#endif /* SWIFT_MINIMAL_MULTI_MAT_HYDRO_H */
diff --git a/src/hydro/MinimalMultiMat/hydro_debug.h b/src/hydro/MinimalMultiMat/hydro_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..d8fe73313fbb175faa9547970c6680346dad0a1b
--- /dev/null
+++ b/src/hydro/MinimalMultiMat/hydro_debug.h
@@ -0,0 +1,53 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2018 Jacob Kegerreis (jacob.kegerreis@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_MINIMAL_MULTI_MAT_HYDRO_DEBUG_H
+#define SWIFT_MINIMAL_MULTI_MAT_HYDRO_DEBUG_H
+
+/**
+ * @file MinimalMultiMat/hydro_debug.h
+ * @brief MinimalMultiMat conservative implementation of SPH (Debugging
+ * routines)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
+ */
+
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+ const struct part* p, const struct xpart* xp) {
+ printf(
+ "x=[%.3e,%.3e,%.3e], "
+ "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
+ "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
+ "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, "
+ "time_bin=%d, mat_id=%d\n",
+ p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
+ xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+ p->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h,
+ p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
+ p->time_bin, p->mat_id);
+}
+
+#endif /* SWIFT_MINIMAL_MULTI_MAT_HYDRO_DEBUG_H */
diff --git a/src/hydro/MinimalMultiMat/hydro_iact.h b/src/hydro/MinimalMultiMat/hydro_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..5984c1c483546d87800792ced0ffcc41e0aaa408
--- /dev/null
+++ b/src/hydro/MinimalMultiMat/hydro_iact.h
@@ -0,0 +1,340 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2018 Jacob Kegerreis (jacob.kegerreis@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_MINIMAL_MULTI_MAT_HYDRO_IACT_H
+#define SWIFT_MINIMAL_MULTI_MAT_HYDRO_IACT_H
+
+/**
+ * @file MinimalMultiMat/hydro_iact.h
+ * @brief MinimalMultiMat conservative implementation of SPH (Neighbour loop
+ * equations)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
+ * Physics, 2012, Volume 231, Issue 3, pp. 759-794.
+ */
+
+#include "adiabatic_index.h"
+#include "const.h"
+#include "minmax.h"
+
+/**
+ * @brief Density interaction between two particles.
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_density(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ float wi, wj, wi_dx, wj_dx;
+
+ /* Get r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ /* Get the masses. */
+ const float mi = pi->mass;
+ const float mj = pj->mass;
+
+ /* Compute density of pi. */
+ const float hi_inv = 1.f / hi;
+ const float ui = r * hi_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ pi->rho += mj * wi;
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+ /* Compute density of pj. */
+ const float hj_inv = 1.f / hj;
+ const float uj = r * hj_inv;
+ kernel_deval(uj, &wj, &wj_dx);
+
+ pj->rho += mi * wj;
+ pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+ pj->density.wcount += wj;
+ pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
+}
+
+/**
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ const struct part *restrict pj, float a, float H) {
+
+ float wi, wi_dx;
+
+ /* Get the masses. */
+ const float mj = pj->mass;
+
+ /* Get r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ const float h_inv = 1.f / hi;
+ const float ui = r * h_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ pi->rho += mj * wi;
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+}
+
+/**
+ * @brief Force interaction between two particles.
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_force(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ /* Cosmological factors entering the EoMs */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
+
+ /* Get r and r inverse. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ /* Recover some data */
+ const float mi = pi->mass;
+ const float mj = pj->mass;
+ const float rhoi = pi->rho;
+ const float rhoj = pj->rho;
+ const float pressurei = pi->force.pressure;
+ const float pressurej = pj->force.pressure;
+
+ /* Get the kernel for hi. */
+ const float hi_inv = 1.0f / hi;
+ const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+ const float xi = r * hi_inv;
+ float wi, wi_dx;
+ kernel_deval(xi, &wi, &wi_dx);
+ const float wi_dr = hid_inv * wi_dx;
+
+ /* Get the kernel for hj. */
+ const float hj_inv = 1.0f / hj;
+ const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
+ const float xj = r * hj_inv;
+ float wj, wj_dx;
+ kernel_deval(xj, &wj, &wj_dx);
+ const float wj_dr = hjd_inv * wj_dx;
+
+ /* Compute gradient terms */
+ const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
+ const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
+
+ /* Compute dv dot r. */
+ const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+ (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
+
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = min(dvdr, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Compute sound speeds and signal velocity */
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
+ const float v_sig = ci + cj - 3.f * mu_ij;
+
+ /* Construct the full viscosity term */
+ const float rho_ij = 0.5f * (rhoi + rhoj);
+ const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+
+ /* Convolve with the kernel */
+ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+
+ /* SPH acceleration term */
+ const float sph_acc_term =
+ (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+
+ /* Assemble the acceleration */
+ const float acc = sph_acc_term + visc_acc_term;
+
+ /* Use the force Luke ! */
+ pi->a_hydro[0] -= mj * acc * dx[0];
+ pi->a_hydro[1] -= mj * acc * dx[1];
+ pi->a_hydro[2] -= mj * acc * dx[2];
+
+ pj->a_hydro[0] += mi * acc * dx[0];
+ pj->a_hydro[1] += mi * acc * dx[1];
+ pj->a_hydro[2] += mi * acc * dx[2];
+
+ /* Get the time derivative for u. */
+ const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
+ const float sph_du_term_j = P_over_rho2_j * dvdr * r_inv * wj_dr;
+
+ /* Viscosity term */
+ const float visc_du_term = 0.5f * visc_acc_term * dvdr;
+
+ /* Assemble the energy equation term */
+ const float du_dt_i = sph_du_term_i + visc_du_term;
+ const float du_dt_j = sph_du_term_j + visc_du_term;
+
+ /* Internal energy time derivatibe */
+ pi->u_dt += du_dt_i * mj;
+ pj->u_dt += du_dt_j * mi;
+
+ /* Get the time derivative for h. */
+ pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+ pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
+
+ /* Update the signal velocity. */
+ pi->force.v_sig = max(pi->force.v_sig, v_sig);
+ pj->force.v_sig = max(pj->force.v_sig, v_sig);
+}
+
+/**
+ * @brief Force interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ const struct part *restrict pj, float a, float H) {
+
+ /* Cosmological factors entering the EoMs */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
+
+ /* Get r and r inverse. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ /* Recover some data */
+ // const float mi = pi->mass;
+ const float mj = pj->mass;
+ const float rhoi = pi->rho;
+ const float rhoj = pj->rho;
+ const float pressurei = pi->force.pressure;
+ const float pressurej = pj->force.pressure;
+
+ /* Get the kernel for hi. */
+ const float hi_inv = 1.0f / hi;
+ const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+ const float xi = r * hi_inv;
+ float wi, wi_dx;
+ kernel_deval(xi, &wi, &wi_dx);
+ const float wi_dr = hid_inv * wi_dx;
+
+ /* Get the kernel for hj. */
+ const float hj_inv = 1.0f / hj;
+ const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
+ const float xj = r * hj_inv;
+ float wj, wj_dx;
+ kernel_deval(xj, &wj, &wj_dx);
+ const float wj_dr = hjd_inv * wj_dx;
+
+ /* Compute gradient terms */
+ const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
+ const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
+
+ /* Compute dv dot r. */
+ const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+ (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
+
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = min(dvdr, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Compute sound speeds and signal velocity */
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
+ const float v_sig = ci + cj - 3.f * mu_ij;
+
+ /* Construct the full viscosity term */
+ const float rho_ij = 0.5f * (rhoi + rhoj);
+ const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+
+ /* Convolve with the kernel */
+ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+
+ /* SPH acceleration term */
+ const float sph_acc_term =
+ (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+
+ /* Assemble the acceleration */
+ const float acc = sph_acc_term + visc_acc_term;
+
+ /* Use the force Luke ! */
+ pi->a_hydro[0] -= mj * acc * dx[0];
+ pi->a_hydro[1] -= mj * acc * dx[1];
+ pi->a_hydro[2] -= mj * acc * dx[2];
+
+ /* Get the time derivative for u. */
+ const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
+
+ /* Viscosity term */
+ const float visc_du_term = 0.5f * visc_acc_term * dvdr;
+
+ /* Assemble the energy equation term */
+ const float du_dt_i = sph_du_term_i + visc_du_term;
+
+ /* Internal energy time derivatibe */
+ pi->u_dt += du_dt_i * mj;
+
+ /* Get the time derivative for h. */
+ pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+
+ /* Update the signal velocity. */
+ pi->force.v_sig = max(pi->force.v_sig, v_sig);
+}
+
+#endif /* SWIFT_MINIMAL_MULTI_MAT_HYDRO_IACT_H */
diff --git a/src/hydro/MinimalMultiMat/hydro_io.h b/src/hydro/MinimalMultiMat/hydro_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..6e55497dc42111de291b8984a7dca048bb774722
--- /dev/null
+++ b/src/hydro/MinimalMultiMat/hydro_io.h
@@ -0,0 +1,200 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2018 Jacob Kegerreis (jacob.kegerreis@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_MINIMAL_MULTI_MAT_HYDRO_IO_H
+#define SWIFT_MINIMAL_MULTI_MAT_HYDRO_IO_H
+
+/**
+ * @file MinimalMultiMat/hydro_io.h
+ * @brief MinimalMultiMat conservative implementation of SPH (i/o routines)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
+ */
+
+#include "adiabatic_index.h"
+#include "hydro.h"
+#include "io_properties.h"
+#include "kernel_hydro.h"
+
+/**
+ * @brief Specifies which particle fields to read from a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+void hydro_read_particles(struct part* parts, struct io_props* list,
+ int* num_fields) {
+
+ *num_fields = 9;
+
+ /* List what we want to read */
+ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, x);
+ list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+ UNIT_CONV_SPEED, parts, v);
+ list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+ parts, mass);
+ list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, h);
+ list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
+ list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, parts, id);
+ list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
+ UNIT_CONV_ACCELERATION, parts, a_hydro);
+ list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
+ UNIT_CONV_DENSITY, parts, rho);
+ list[8] =
+ io_make_input_field("MaterialID", INT, 1, OPTIONAL, 1, parts, mat_id);
+}
+
+void convert_P(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ ret[0] = hydro_get_comoving_pressure(p);
+}
+
+void convert_part_pos(const struct engine* e, const struct part* p,
+ const struct xpart* xp, double* ret) {
+
+ if (e->s->periodic) {
+ ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
+ ret[1] = box_wrap(p->x[1], 0.0, e->s->dim[1]);
+ ret[2] = box_wrap(p->x[2], 0.0, e->s->dim[2]);
+ } else {
+ ret[0] = p->x[0];
+ ret[1] = p->x[1];
+ ret[2] = p->x[2];
+ }
+}
+
+void convert_part_vel(const struct engine* e, const struct part* p,
+ const struct xpart* xp, 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, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ 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);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a2_inv;
+ ret[1] *= cosmo->a2_inv;
+ ret[2] *= cosmo->a2_inv;
+}
+
+void convert_part_potential(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ if (p->gpart != NULL)
+ ret[0] = gravity_get_comoving_potential(p->gpart);
+ else
+ ret[0] = 0.f;
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param xparts The extended particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+ struct io_props* list, int* num_fields) {
+
+ *num_fields = 10;
+
+ /* List what we want to write */
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
+ list[2] =
+ io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
+ list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+ parts, h);
+ list[4] = io_make_output_field("InternalEnergy", FLOAT, 1,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
+ list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
+ UNIT_CONV_NO_UNITS, parts, id);
+ list[6] =
+ io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
+ list[7] = io_make_output_field("MaterialID", INT, 1, UNIT_CONV_NO_UNITS,
+ parts, mat_id);
+ list[8] = io_make_output_field_convert_part(
+ "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, xparts, convert_P);
+ list[9] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+ UNIT_CONV_POTENTIAL, parts,
+ xparts, convert_part_potential);
+}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+void hydro_write_flavour(hid_t h_grpsph) {
+
+ /* Viscosity and thermal conduction */
+ /* Nothing in this minimal model... */
+ io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
+ io_write_attribute_s(h_grpsph, "Viscosity Model",
+ "Minimal treatment as in Monaghan (1992)");
+
+ /* Time integration properties */
+ io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
+ const_max_u_change);
+}
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+int writeEntropyFlag() { return 0; }
+
+#endif /* SWIFT_MINIMAL_MULTI_MAT_HYDRO_IO_H */
diff --git a/src/hydro/MinimalMultiMat/hydro_part.h b/src/hydro/MinimalMultiMat/hydro_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..dad13e889aa531636e34846825109086177b3dae
--- /dev/null
+++ b/src/hydro/MinimalMultiMat/hydro_part.h
@@ -0,0 +1,180 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2018 Jacob Kegerreis (jacob.kegerreis@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_MINIMAL_MULTI_MAT_HYDRO_PART_H
+#define SWIFT_MINIMAL_MULTI_MAT_HYDRO_PART_H
+
+/**
+ * @file MinimalMultiMat/hydro_part.h
+ * @brief MinimalMultiMat conservative implementation of SPH (Particle
+ * definition)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
+ * Physics, 2012, Volume 231, Issue 3, pp. 759-794.
+ */
+
+#include "chemistry_struct.h"
+#include "cooling_struct.h"
+#include "equation_of_state.h" // For enum material_id
+
+/**
+ * @brief Particle fields not needed during the SPH loops over neighbours.
+ *
+ * This structure contains the particle fields that are not used in the
+ * density or force loops. Quantities should be used in the kick, drift and
+ * potentially ghost tasks only.
+ */
+struct xpart {
+
+ /*! Offset between current position and position at last tree rebuild. */
+ float x_diff[3];
+
+ /*! Offset between the current position and position at the last sort. */
+ float x_diff_sort[3];
+
+ /*! Velocity at the last full step. */
+ float v_full[3];
+
+ /*! Gravitational acceleration at the last full step. */
+ float a_grav[3];
+
+ /*! Internal energy at the last full step. */
+ float u_full;
+
+ /*! Additional data used to record cooling information */
+ struct cooling_xpart_data cooling_data;
+
+} SWIFT_STRUCT_ALIGN;
+
+/**
+ * @brief Particle fields for the SPH particles
+ *
+ * The density and force substructures are used to contain variables only used
+ * within the density and force loops over neighbours. All more permanent
+ * variables should be declared in the main part of the part structure,
+ */
+struct part {
+
+ /*! Particle unique ID. */
+ long long id;
+
+ /*! Pointer to corresponding gravity part. */
+ struct gpart* gpart;
+
+ /*! Particle position. */
+ double x[3];
+
+ /*! Particle predicted velocity. */
+ float v[3];
+
+ /*! Particle acceleration. */
+ float a_hydro[3];
+
+ /*! Particle mass. */
+ float mass;
+
+ /*! Particle smoothing length. */
+ float h;
+
+ /*! Particle internal energy. */
+ float u;
+
+ /*! Time derivative of the internal energy. */
+ float u_dt;
+
+ /*! Particle density. */
+ float rho;
+
+ /* Store density/force specific stuff. */
+ union {
+
+ /**
+ * @brief Structure for the variables only used in the density loop over
+ * neighbours.
+ *
+ * Quantities in this sub-structure should only be accessed in the density
+ * loop over neighbours and the ghost task.
+ */
+ struct {
+
+ /*! Neighbour number count. */
+ float wcount;
+
+ /*! Derivative of the neighbour number with respect to h. */
+ float wcount_dh;
+
+ /*! Derivative of density with respect to h */
+ float rho_dh;
+
+ } density;
+
+ /**
+ * @brief Structure for the variables only used in the force loop over
+ * neighbours.
+ *
+ * Quantities in this sub-structure should only be accessed in the force
+ * loop over neighbours and the ghost, drift and kick tasks.
+ */
+ struct {
+
+ /*! "Grad h" term */
+ float f;
+
+ /*! Particle pressure. */
+ float pressure;
+
+ /*! Particle soundspeed. */
+ float soundspeed;
+
+ /*! Particle signal velocity */
+ float v_sig;
+
+ /*! Time derivative of smoothing length */
+ float h_dt;
+
+ } force;
+ };
+
+ /* Chemistry information */
+ struct chemistry_part_data chemistry_data;
+
+ /*! Material identifier flag */
+ enum eos_planetary_material_id mat_id;
+
+ /*! Time-step length */
+ timebin_t time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Time of the last drift */
+ integertime_t ti_drift;
+
+ /* Time of the last kick */
+ integertime_t ti_kick;
+
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_MINIMAL_MULTI_MAT_HYDRO_PART_H */
diff --git a/src/hydro_io.h b/src/hydro_io.h
index 639c2f3ae640d7b74e6a2507bd4e3d5ad5625171..c5f9aae9f7bb5f581b09dcd3f309c5fa95f33e51 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -35,6 +35,8 @@
#include "./hydro/Gizmo/hydro_io.h"
#elif defined(SHADOWFAX_SPH)
#include "./hydro/Shadowswift/hydro_io.h"
+#elif defined(MINIMAL_MULTI_MAT_SPH)
+#include "./hydro/MinimalMultiMat/hydro_io.h"
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/part.h b/src/part.h
index 3e9ab5180e7cbfb4622b4a05d396a7f785f08b0b..0cfe731e52167825b89c3e732620bf47ba75bdad 100644
--- a/src/part.h
+++ b/src/part.h
@@ -62,6 +62,9 @@
#include "./hydro/Shadowswift/hydro_part.h"
#define hydro_need_extra_init_loop 0
#define EXTRA_HYDRO_LOOP
+#elif defined(MINIMAL_MULTI_MAT_SPH)
+#include "./hydro/MinimalMultiMat/hydro_part.h"
+#define hydro_need_extra_init_loop 0
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/runner.c b/src/runner.c
index c275225e3a4bb881545571110795a5beb17bf6e8..5a382a916abda9dddd45dcb5e703578457f0548a 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -847,8 +847,9 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
error("Smoothing length failed to converge on %i particles.", count);
}
#else
- if (count)
+ if (count) {
error("Smoothing length failed to converge on %i particles.", count);
+ }
#endif
/* Be clean */
diff --git a/tests/Makefile.am b/tests/Makefile.am
index ca80b30254ac82b7fa26220d706713974bcc1a20..fb3384fa56936b3c7c1fe1b415a2201e419ba566 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -1,16 +1,16 @@
# This file is part of SWIFT.
# Copyright (c) 2015 matthieu.schaller@durham.ac.uk.
-#
+#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU 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 General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
@@ -27,7 +27,7 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
testMatrixInversion testThreadpool testDump testLogger testInteractions.sh \
testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
- testPotentialPair
+ testPotentialPair testEOS
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
@@ -37,7 +37,7 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testAdiabaticIndex testRiemannExact testRiemannTRRS \
testRiemannHLLC testMatrixInversion testDump testLogger \
testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
- testGravityDerivatives testPotentialSelf testPotentialPair
+ testGravityDerivatives testPotentialSelf testPotentialPair testEOS
# Rebuild tests when SWIFT is updated.
$(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -105,6 +105,8 @@ testPotentialSelf_SOURCES = testPotentialSelf.c
testPotentialPair_SOURCES = testPotentialPair.c
+testEOS_SOURCES = testEOS.c
+
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testActivePair.sh \
test27cells.sh test27cellsPerturbed.sh testParser.sh testPeriodicBC.sh \
@@ -112,4 +114,5 @@ EXTRA_DIST = testReading.sh makeInput.py testActivePair.sh \
difffloat.py tolerance_125_normal.dat tolerance_125_perturbed.dat \
tolerance_27_normal.dat tolerance_27_perturbed.dat tolerance_27_perturbed_h.dat tolerance_27_perturbed_h2.dat \
tolerance_testInteractions.dat tolerance_pair_active.dat tolerance_pair_force_active.dat \
- fft_params.yml tolerance_periodic_BC_normal.dat tolerance_periodic_BC_perturbed.dat
+ fft_params.yml tolerance_periodic_BC_normal.dat tolerance_periodic_BC_perturbed.dat \
+ testEOS.sh testEOS_plot.sh
diff --git a/tests/test125cells.c b/tests/test125cells.c
index ddebce8cabda227e96e04996214104354703c185..1dfd5ea73a1a40c1106999addf7f7ef6352b79d9 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -120,6 +120,8 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(MINIMAL_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
+#elif defined(MINIMAL_MULTI_MAT_SPH)
+ part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
part->primitives.P = pressure;
#else
@@ -403,7 +405,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
main_cell->parts[pid].v[2], main_cell->parts[pid].h,
hydro_get_comoving_density(&main_cell->parts[pid]),
-#if defined(MINIMAL_SPH) || defined(SHADOWFAX_SPH)
+#if defined(MINIMAL_SPH) || defined(MINIMAL_MULTI_MAT_SPH) || \
+ defined(SHADOWFAX_SPH)
0.f,
#else
main_cell->parts[pid].density.div_v,
diff --git a/tests/testEOS.c b/tests/testEOS.c
new file mode 100644
index 0000000000000000000000000000000000000000..2e72d3d1768f3a6ea4ab1665a099efeb28f8f3f9
--- /dev/null
+++ b/tests/testEOS.c
@@ -0,0 +1,278 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2018 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <fenv.h>
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+#include <pthread.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+
+/* Conditional headers. */
+#ifdef HAVE_LIBZ
+#include <zlib.h>
+#endif
+
+/* Local headers. */
+#include "equation_of_state.h"
+#include "swift.h"
+
+/* Engine policy flags. */
+#ifndef ENGINE_POLICY
+#define ENGINE_POLICY engine_policy_none
+#endif
+
+/**
+ * @brief Write a list of densities, energies, and resulting pressures to file
+ * from an equation of state.
+ *
+ * WORK IN PROGRESS
+ *
+ * So far only does the Hubbard & MacFarlane (1980) equations of state.
+ *
+ * Usage:
+ * $ ./testEOS (mat_id) (do_output)
+ *
+ * Sys args (optional):
+ * mat_id | int | Material ID, see equation_of_state.h for the options.
+ * Default: 201 (= id_HM80_ice).
+ *
+ * do_output | int | Set 1 to write the output file of rho, u, P values,
+ * set 0 for no output. Default: 0.
+ *
+ * Output text file contains:
+ * header
+ * num_rho num_u mat_id # Header values
+ * rho_0 rho_1 rho_2 ... rho_num_rho # Array of densities, rho
+ * u_0 u_1 u_2 ... u_num_u # Array of energies, u
+ * P_0_0 P_0_1 ... P_0_num_u # Array of pressures, P(rho, u)
+ * P_1_0 ... ... P_1_num_u
+ * ... ... ... ...
+ * P_num_rho_0 ... P_num_rho_num_u
+ *
+ * Note that the values tested extend beyond the range that most EOS are
+ * designed for (e.g. outside table limits), to help test the EOS in case of
+ * unexpected particle behaviour.
+ *
+ */
+
+#ifdef EOS_PLANETARY
+int main(int argc, char *argv[]) {
+ float rho, log_rho, log_u, P;
+ struct unit_system us;
+ const struct phys_const *phys_const = 0; // Unused placeholder
+ const struct swift_params *params = 0; // Unused placeholder
+ const float J_kg_to_erg_g = 1e4; // Convert J/kg to erg/g
+ char filename[64];
+ // Output table params
+ const int num_rho = 100, num_u = 100;
+ float log_rho_min = logf(1e-4), log_rho_max = logf(30.f),
+ log_u_min = logf(1e4), log_u_max = logf(1e10),
+ log_rho_step = (log_rho_max - log_rho_min) / (num_rho - 1.f),
+ log_u_step = (log_u_max - log_u_min) / (num_u - 1.f);
+ float A1_rho[num_rho], A1_u[num_u];
+ // Sys args
+ int mat_id, do_output;
+ // Default sys args
+ const int mat_id_def = eos_planetary_id_HM80_ice;
+ const int do_output_def = 0;
+
+ // Check the number of system arguments and use defaults if not provided
+ switch (argc) {
+ case 1:
+ // Default both
+ mat_id = mat_id_def;
+ do_output = do_output_def;
+ break;
+
+ case 2:
+ // Read mat_id, default do_output
+ mat_id = atoi(argv[1]);
+ do_output = do_output_def;
+ break;
+
+ case 3:
+ // Read both
+ mat_id = atoi(argv[1]);
+ do_output = atoi(argv[2]);
+ break;
+
+ default:
+ error("Invalid number of system arguments!\n");
+ mat_id = mat_id_def; // Ignored, just here to keep the compiler happy
+ do_output = do_output_def;
+ };
+
+ /* Greeting message */
+ printf("This is %s\n", package_description());
+
+ // Check material ID
+ // Material base type
+ switch ((int)(mat_id / eos_planetary_type_factor)) {
+ // Tillotson
+ case eos_planetary_type_Til:
+ switch (mat_id) {
+ case eos_planetary_id_Til_iron:
+ printf(" Tillotson iron \n");
+ break;
+
+ case eos_planetary_id_Til_granite:
+ printf(" Tillotson granite \n");
+ break;
+
+ case eos_planetary_id_Til_water:
+ printf(" Tillotson water \n");
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d \n", mat_id);
+ };
+ break;
+
+ // Hubbard & MacFarlane (1980)
+ case eos_planetary_type_HM80:
+ switch (mat_id) {
+ case eos_planetary_id_HM80_HHe:
+ printf(" Hubbard & MacFarlane (1980) hydrogen-helium atmosphere \n");
+ break;
+
+ case eos_planetary_id_HM80_ice:
+ printf(" Hubbard & MacFarlane (1980) ice mix \n");
+ break;
+
+ case eos_planetary_id_HM80_rock:
+ printf(" Hubbard & MacFarlane (1980) rock mix \n");
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d \n", mat_id);
+ };
+ break;
+
+ // ANEOS
+ case eos_planetary_type_ANEOS:
+ switch (mat_id) {
+ case eos_planetary_id_ANEOS_iron:
+ printf(" ANEOS iron \n");
+ break;
+
+ case eos_planetary_id_MANEOS_forsterite:
+ printf(" MANEOS forsterite \n");
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d \n", mat_id);
+ };
+ break;
+
+ // SESAME
+ case eos_planetary_type_SESAME:
+ switch (mat_id) {
+ case eos_planetary_id_SESAME_iron:
+ printf(" SESAME iron \n");
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d \n", mat_id);
+ };
+ break;
+
+ default:
+ error("Unknown material type! mat_id = %d \n", mat_id);
+ }
+
+ // Convert to internal units (Earth masses and radii)
+ units_init(&us, 5.9724e27, 6.3710e8, 1.f, 1.f, 1.f);
+ log_rho_min -= logf(units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY));
+ log_rho_max -= logf(units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY));
+ log_u_min += logf(J_kg_to_erg_g / units_cgs_conversion_factor(
+ &us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
+ log_u_max += logf(J_kg_to_erg_g / units_cgs_conversion_factor(
+ &us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
+
+ // Initialise the EOS materials
+ eos_init(&eos, phys_const, &us, params);
+
+ // Output file
+ sprintf(filename, "testEOS_rho_u_P_%d.txt", mat_id);
+ FILE *f = fopen(filename, "w");
+ if (f == NULL) {
+ printf("Could not open output file!\n");
+ exit(EXIT_FAILURE);
+ }
+
+ if (do_output == 1) {
+ fprintf(f, "Density Sp.Int.Energy mat_id \n");
+ fprintf(f, "%d %d %d \n", num_rho, num_u, mat_id);
+ }
+
+ // Densities
+ log_rho = log_rho_min;
+ for (int i = 0; i < num_rho; i++) {
+ A1_rho[i] = exp(log_rho);
+ log_rho += log_rho_step;
+
+ if (do_output == 1)
+ fprintf(f, "%.6e ",
+ A1_rho[i] * units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY));
+ }
+ if (do_output == 1) fprintf(f, "\n");
+
+ // Sp. int. energies
+ log_u = log_u_min;
+ for (int i = 0; i < num_u; i++) {
+ A1_u[i] = exp(log_u);
+ log_u += log_u_step;
+
+ if (do_output == 1)
+ fprintf(f, "%.6e ", A1_u[i] * units_cgs_conversion_factor(
+ &us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
+ }
+ if (do_output == 1) fprintf(f, "\n");
+
+ // Pressures
+ for (int i = 0; i < num_rho; i++) {
+ rho = A1_rho[i];
+
+ for (int j = 0; j < num_u; j++) {
+ P = gas_pressure_from_internal_energy(rho, A1_u[j], mat_id);
+
+ if (do_output == 1)
+ fprintf(f, "%.6e ",
+ P * units_cgs_conversion_factor(&us, UNIT_CONV_PRESSURE));
+ }
+
+ if (do_output == 1) fprintf(f, "\n");
+ }
+ fclose(f);
+
+ return 0;
+}
+#else
+int main() { return 0; }
+#endif
diff --git a/tests/testEOS.py b/tests/testEOS.py
new file mode 100644
index 0000000000000000000000000000000000000000..363bab200b58c65fa24cc033af4b8d3c04b7b503
--- /dev/null
+++ b/tests/testEOS.py
@@ -0,0 +1,176 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 2018 Jacob Kegerreis (jacob.kegerreis@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/>.
+ #
+ ##############################################################################
+"""
+Plot the output of testEOS to show how the equation of state pressure varies
+with density and specific internal energy.
+
+Usage:
+ python testEOS.py (mat_id)
+
+ Sys args (optional):
+ mat_id | int | Material ID, see equation_of_state.h for the options.
+ Default: 201 (= HM80_ice).
+
+Text file contains:
+ header
+ num_rho num_u mat_id # Header info
+ rho_0 rho_1 rho_2 ... rho_num_rho # Array of densities, rho
+ u_0 u_1 u_2 ... u_num_u # Array of energies, u
+ P_0_0 P_0_1 ... P_0_num_u # Array of pressures, P(rho, u)
+ P_1_0 ... ... P_1_num_u
+ ... ... ... ...
+ P_num_rho_0 ... P_num_rho_num_u
+
+Note that the values tested extend beyond the range that most EOS are
+designed for (e.g. outside table limits), to help test the EOS in case of
+unexpected particle behaviour.
+"""
+
+# ========
+# Modules and constants
+# ========
+from __future__ import division
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import sys
+
+# Material types (copied from src/equation_of_state/planetary/equation_of_state.h)
+type_factor = 100
+Di_type = {
+ 'Til' : 1,
+ 'HM80' : 2,
+ 'ANEOS' : 3,
+ 'SESAME' : 4,
+}
+Di_material = {
+ # Tillotson
+ 'Til_iron' : Di_type['Til']*type_factor,
+ 'Til_granite' : Di_type['Til']*type_factor + 1,
+ 'Til_water' : Di_type['Til']*type_factor + 2,
+ # Hubbard & MacFarlane (1980) Uranus/Neptune
+ 'HM80_HHe' : Di_type['HM80']*type_factor, # Hydrogen-helium atmosphere
+ 'HM80_ice' : Di_type['HM80']*type_factor + 1, # H20-CH4-NH3 ice mix
+ 'HM80_rock' : Di_type['HM80']*type_factor + 2, # SiO2-MgO-FeS-FeO rock mix
+ # ANEOS
+ 'ANEOS_iron' : Di_type['ANEOS']*type_factor,
+ 'MANEOS_forsterite' : Di_type['ANEOS']*type_factor + 1,
+ # SESAME
+ 'SESAME_iron' : Di_type['SESAME']*type_factor,
+}
+# Invert so the mat_id are the keys
+Di_mat_id = {mat_id : mat for mat, mat_id in Di_material.iteritems()}
+
+# Unit conversion
+Ba_to_Mbar = 1e-12
+erg_g_to_J_kg = 1e-4
+
+if __name__ == '__main__':
+ # Sys args
+ try:
+ mat_id = int(sys.argv[1])
+ except IndexError:
+ mat_id = Di_material['HM80_ice']
+
+ # Check the material
+ try:
+ mat = Di_mat_id[mat_id]
+ print mat
+ sys.stdout.flush()
+ except KeyError:
+ print "Error: unknown material ID! mat_id = %d" % mat_id
+ print "Materials:"
+ for mat_id, mat in sorted(Di_mat_id.iteritems()):
+ print " %s%s%d" % (mat, (20 - len("%s" % mat))*" ", mat_id)
+
+ filename = "testEOS_rho_u_P_%d.txt" % mat_id
+
+ # Load the header info and density and energy arrays
+ with open(filename) as f:
+ f.readline()
+ num_rho, num_u, mat_id = np.array(f.readline().split(), dtype=int)
+ A1_rho = np.array(f.readline().split(), dtype=float)
+ A1_u = np.array(f.readline().split(), dtype=float)
+
+ # Load pressure array
+ A2_P = np.loadtxt(filename, skiprows=4)
+
+ # Convert pressures from cgs Barye to Mbar
+ A2_P *= Ba_to_Mbar
+ # Convert energies from cgs to SI
+ A1_u *= erg_g_to_J_kg
+
+ # Check that the numbers are right
+ assert A1_rho.shape == (num_rho,)
+ assert A1_u.shape == (num_u,)
+ assert A2_P.shape == (num_rho, num_u)
+
+ # Plot
+ plt.figure(figsize=(7, 7))
+ ax = plt.gca()
+
+ # P(rho) at fixed u
+ num_u_fix = 9
+ A1_idx = np.floor(np.linspace(0, num_u - 1, num_u_fix)).astype(int)
+ A1_colour = matplotlib.cm.rainbow(np.linspace(0, 1, num_u_fix))
+
+ for i, idx in enumerate(A1_idx):
+ plt.plot(A1_rho, A2_P[:, idx], c=A1_colour[i],
+ label=r"%.2e" % A1_u[idx])
+
+ plt.legend(title="Sp. Int. Energy (J kg$^{-1}$)")
+ plt.xscale('log')
+ plt.yscale('log')
+ plt.xlabel(r"Density (g cm$^{-3}$)")
+ plt.ylabel(r"Pressure (Mbar)")
+ plt.title(mat)
+ plt.tight_layout()
+
+ plt.savefig("testEOS_%d.png" % mat_id)
+ plt.close()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testEOS.sh b/tests/testEOS.sh
new file mode 100755
index 0000000000000000000000000000000000000000..411ac746be186bfe5758e03c2a852e081daefd10
--- /dev/null
+++ b/tests/testEOS.sh
@@ -0,0 +1,23 @@
+#!/bin/bash
+
+echo ""
+
+rm -f testEOS_rho_u_P_*.txt
+
+echo "Running testEOS for each planetary material"
+
+A1_mat_id=(
+ 100
+ 101
+ 102
+ 200
+ 201
+ 202
+)
+
+for mat_id in "${A1_mat_id[@]}"
+do
+ ./testEOS "$mat_id" 1
+done
+
+exit $?
diff --git a/tests/testEOS_plot.sh b/tests/testEOS_plot.sh
new file mode 100755
index 0000000000000000000000000000000000000000..39108c5e19d8f4474de508e205951a1fd0aebcc9
--- /dev/null
+++ b/tests/testEOS_plot.sh
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+echo ""
+
+echo "Plotting testEOS output for each planetary material"
+
+A1_mat_id=(
+ 100
+ 101
+ 102
+ 200
+ 201
+ 202
+)
+
+for mat_id in "${A1_mat_id[@]}"
+do
+ python ./testEOS.py "$mat_id"
+done
+
+exit $?
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index 9c5dd36415970ff2a53220aa56cecba6fc5fe193..e1581c52863072b9460e764eb65597fff5817a49 100644
--- a/tests/testInteractions.c
+++ b/tests/testInteractions.c
@@ -104,7 +104,8 @@ struct part *make_particles(size_t count, double *offset, double spacing,
*/
void prepare_force(struct part *parts, size_t count) {
-#if !defined(GIZMO_SPH) && !defined(SHADOWFAX_SPH) && !defined(MINIMAL_SPH)
+#if !defined(GIZMO_SPH) && !defined(SHADOWFAX_SPH) && !defined(MINIMAL_SPH) && \
+ !defined(MINIMAL_MULTI_MAT_SPH)
struct part *p;
for (size_t i = 0; i < count; ++i) {
p = &parts[i];
@@ -131,7 +132,8 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
"%8.5f %8.5f %13e %13e %13e %13e %13e %8.5f %8.5f\n",
p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->h,
hydro_get_comoving_density(p),
-#if defined(MINIMAL_SPH) || defined(SHADOWFAX_SPH)
+#if defined(MINIMAL_SPH) || defined(MINIMAL_MULTI_MAT_SPH) || \
+ defined(SHADOWFAX_SPH)
0.f,
#else
p->density.div_v,