diff --git a/.gitignore b/.gitignore
index 27ff2cf7ff880de2007e27978fcc1cfa22d2bd75..d2f69ad4e8759839082a38d4558a5ab7bac0e18d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -165,6 +165,8 @@ m4/lt~obsolete.m4
/stamp-h1
/test-driver
+src/equation_of_state/planetary/*.txt
+
# Intel compiler optimization reports
*.optrpt
@@ -316,3 +318,6 @@ sympy-plots-for-*.tex/
#ctags
*tags
+
+# vim
+*.swp
diff --git a/configure.ac b/configure.ac
index 56f88116c845def33300b87630524f06fec666bb..6ec1ece5ed2fef20ceab484eff4e09fce808e635 100644
--- a/configure.ac
+++ b/configure.ac
@@ -576,7 +576,7 @@ if test "x$with_profiler" != "xno"; then
proflibs="-lprofiler"
fi
AC_CHECK_LIB([profiler],[ProfilerFlush],
- [have_profiler="yes"
+ [have_profiler="yes"
AC_DEFINE([WITH_PROFILER],1,[Link against the gperftools profiling library.])],
[have_profiler="no"], $proflibs)
@@ -973,7 +973,7 @@ esac
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
- [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, default, gizmo-mfv, gizmo-mfm, shadowfax, minimal-multi-mat, debug default: gadget2@:>@]
+ [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, debug default: gadget2@:>@]
)],
[with_hydro="$withval"],
[with_hydro="gadget2"]
@@ -1012,10 +1012,11 @@ 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])
+ planetary)
+ AC_DEFINE([PLANETARY_SPH], [1], [Planetary SPH])
;;
+
*)
AC_MSG_ERROR([Unknown hydrodynamics scheme: $with_hydro])
;;
diff --git a/examples/MoonFormingImpact/README.md b/examples/MoonFormingImpact/README.md
deleted file mode 100644
index 97a84f67c6aeeff4176a1385381f1cfe9e340c91..0000000000000000000000000000000000000000
--- a/examples/MoonFormingImpact/README.md
+++ /dev/null
@@ -1,34 +0,0 @@
-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
deleted file mode 100755
index 7d63943c2c5dc3bd4ab88e63a2abba62cc3f04a5..0000000000000000000000000000000000000000
--- a/examples/MoonFormingImpact/get_init_cond.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-#!/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
deleted file mode 100644
index 323adf7f3ac73f41b45b50eaa76a95033dca35d7..0000000000000000000000000000000000000000
--- a/examples/MoonFormingImpact/moon_forming_impact.yml
+++ /dev/null
@@ -1,48 +0,0 @@
-# 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
-
-# Parameters related to the equation of state
-EoS:
- planetary_use_Til: 1 # Whether to prepare the Tillotson EOS
diff --git a/examples/MoonFormingImpact/plot.py b/examples/MoonFormingImpact/plot.py
deleted file mode 100644
index aa0d64a5d0d06709d51b1db231c507e22861f36c..0000000000000000000000000000000000000000
--- a/examples/MoonFormingImpact/plot.py
+++ /dev/null
@@ -1,285 +0,0 @@
-"""
-###############################################################################
-# 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
deleted file mode 100755
index 165dae3a24a9c30960959fbb37aa6e1da2eb851f..0000000000000000000000000000000000000000
--- a/examples/MoonFormingImpact/run.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-#!/bin/bash
-../swift -G -s -t 8 moon_forming_impact.yml
diff --git a/examples/UranusImpact/README.md b/examples/UranusImpact/README.md
deleted file mode 100644
index 178a3937ecbe527df8e8e82a0d8fd8bcbf9dbef7..0000000000000000000000000000000000000000
--- a/examples/UranusImpact/README.md
+++ /dev/null
@@ -1,40 +0,0 @@
-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
deleted file mode 100755
index e12e009adfbd727cb2452ac21c477b3ecd77b9c9..0000000000000000000000000000000000000000
--- a/examples/UranusImpact/get_init_cond.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-#!/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
deleted file mode 100644
index 3db3bf21bb15862ec524a069c38e47564b48df1d..0000000000000000000000000000000000000000
--- a/examples/UranusImpact/plot.py
+++ /dev/null
@@ -1,291 +0,0 @@
-"""
-###############################################################################
-# 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
deleted file mode 100755
index c6773b7e40fff3fa312dfcb5ba4ada9d9e4b1b8d..0000000000000000000000000000000000000000
--- a/examples/UranusImpact/run.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-#!/bin/bash
-../swift -G -s -t 8 uranus_impact.yml
diff --git a/examples/UranusImpact/uranus_impact.yml b/examples/UranusImpact/uranus_impact.yml
deleted file mode 100644
index fabddca00f80fcdd79ff6114ff0544cd251046f4..0000000000000000000000000000000000000000
--- a/examples/UranusImpact/uranus_impact.yml
+++ /dev/null
@@ -1,51 +0,0 @@
-# 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
-
-# Parameters related to the equation of state
-EoS:
- planetary_use_HM80: 1 # Whether to prepare the Hubbard & MacFarlane (1980) EOS
- # Table file paths
- planetary_HM80_HHe_table_file: /gpfs/data/dc-kege1/gihr_data/P_rho_u_HHe.txt
- planetary_HM80_ice_table_file: /gpfs/data/dc-kege1/gihr_data/P_rho_u_ice.txt
- planetary_HM80_rock_table_file: /gpfs/data/dc-kege1/gihr_data/P_rho_u_roc.txt
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index ddb71c594122a3e8d6ddbd7c5b73e0474b404a75..5fb48eb17d1c210d2a320917e1fd5d1ad67ddd94 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -75,6 +75,7 @@ Snapshots:
time_first: 0. # (Optional) Time of the first output if non-cosmological time-integration (in internal units)
delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
compression: 0 # (Optional) Set the level of compression of the HDF5 datasets [0-9]. 0 does no compression.
+ label_first: 0 # (Optional) An additional offset for the snapshot output label
label_delta: 1 # (Optional) Set the integer increment between snapshot output labels
UnitMass_in_cgs: 1 # (Optional) Unit system for the outputs (Grams)
UnitLength_in_cgs: 1 # (Optional) Unit system for the outputs (Centimeters)
@@ -142,9 +143,13 @@ EoS:
planetary_use_ANEOS: 0 # (Optional) Whether to prepare the ANEOS EOS
planetary_use_SESAME: 0 # (Optional) Whether to prepare the SESAME EOS
# (Optional) Table file paths
- planetary_HM80_HHe_table_file: HM80_HHe.txt
- planetary_HM80_ice_table_file: HM80_ice.txt
- planetary_HM80_rock_table_file: HM80_rock.txt
+ planetary_HM80_HHe_table_file: ./equation_of_state/planetary_HM80_HHe.txt
+ planetary_HM80_ice_table_file: ./equation_of_state/planetary_HM80_ice.txt
+ planetary_HM80_rock_table_file: ./equation_of_state/planetary_HM80_rock.txt
+ planetary_SESAME_iron_table_file: ./equation_of_state/planetary_SESAME_iron_2140.txt
+ planetary_SESAME_basalt_table_file: ./equation_of_state/planetary_SESAME_basalt_7530.txt
+ planetary_SESAME_water_table_file: ./equation_of_state/planetary_SESAME_water_7154.txt
+ planetary_SS08_water_table_file: ./equation_of_state/planetary_SS08_water.txt
# Parameters related to external potentials --------------------------------------------
diff --git a/src/cell.c b/src/cell.c
index 85f8531c261ed8878bff4e32ba2419616b754372..05a1990dd5355877cadcd0526252a6266dddcf6b 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -62,6 +62,7 @@
#include "space.h"
#include "space_getsid.h"
#include "timers.h"
+#include "tools.h"
/* Global variables. */
int cell_next_tag = 0;
@@ -2551,6 +2552,36 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
}
#endif
+#ifdef PLANETARY_SPH
+ /* Remove particles that cross the non-periodic box edge */
+ if (!(e->s->periodic)) {
+ for (int i = 0; i < 3; i++) {
+ if ((p->x[i] - xp->v_full[i] * dt_drift > e->s->dim[i]) ||
+ (p->x[i] - xp->v_full[i] * dt_drift < 0.f) ||
+ ((p->mass != 0.f) && ((p->x[i] < 0.01f * e->s->dim[i]) ||
+ (p->x[i] > 0.99f * e->s->dim[i])))) {
+ /* (TEMPORARY) Crudely stop the particle manually */
+ message(
+ "Particle %lld hit a box edge. \n"
+ " pos=%.4e %.4e %.4e vel=%.2e %.2e %.2e",
+ p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2]);
+ for (int j = 0; j < 3; j++) {
+ p->v[j] = 0.f;
+ p->gpart->v_full[j] = 0.f;
+ xp->v_full[j] = 0.f;
+ }
+ p->h = hydro_h_max;
+ p->time_bin = time_bin_inhibited;
+ p->gpart->time_bin = time_bin_inhibited;
+ hydro_part_has_no_neighbours(p, xp, e->cosmology);
+ p->mass = 0.f;
+ p->gpart->mass = 0.f;
+ break;
+ }
+ }
+ }
+#endif
+
/* Limit h to within the allowed range */
p->h = min(p->h, hydro_h_max);
@@ -2653,6 +2684,26 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
/* Drift... */
drift_gpart(gp, dt_drift, ti_old_gpart, ti_current);
+#ifdef PLANETARY_SPH
+ /* Remove particles that cross the non-periodic box edge */
+ if (!(e->s->periodic)) {
+ for (int i = 0; i < 3; i++) {
+ if ((gp->x[i] - gp->v_full[i] * dt_drift > e->s->dim[i]) ||
+ (gp->x[i] - gp->v_full[i] * dt_drift < 0.f) ||
+ ((gp->mass != 0.f) && ((gp->x[i] < 0.01f * e->s->dim[i]) ||
+ (gp->x[i] > 0.99f * e->s->dim[i])))) {
+ /* (TEMPORARY) Crudely stop the particle manually */
+ for (int j = 0; j < 3; j++) {
+ gp->v_full[j] = 0.f;
+ }
+ gp->time_bin = time_bin_inhibited;
+ gp->mass = 0.f;
+ break;
+ }
+ }
+ }
+#endif
+
/* Init gravity force fields. */
if (gpart_is_active(gp, e)) {
gravity_init_gpart(gp);
diff --git a/src/debug.c b/src/debug.c
index da8ef0e118b57a6aa94577898b03bcf7c56b006a..5c2523da7b97e44febe52405e954eef1c5f865d2 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -56,8 +56,8 @@
#include "./hydro/GizmoMFM/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"
+#elif defined(PLANETARY_SPH)
+#include "./hydro/Planetary/hydro_debug.h"
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/engine.c b/src/engine.c
index 8ba089295da4b52d8d094c7740faeb701bb61783..d0668fa3c33f2a17ae112216d3b33d031376277e 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -62,6 +62,7 @@
#include "cosmology.h"
#include "cycle.h"
#include "debug.h"
+#include "equation_of_state.h"
#include "error.h"
#include "gravity.h"
#include "gravity_cache.h"
@@ -5767,6 +5768,10 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
parser_get_param_string(params, "Snapshots:basename", e->snapshot_base_name);
e->snapshot_compression =
parser_get_opt_param_int(params, "Snapshots:compression", 0);
+ e->snapshot_label_first =
+ parser_get_opt_param_int(params, "Snapshots:label_first", 0);
+ if (e->snapshot_label_first < 0)
+ error("Snapshots:label_first must be zero or positive");
e->snapshot_label_delta =
parser_get_opt_param_int(params, "Snapshots:label_delta", 1);
e->snapshot_units = (struct unit_system *)malloc(sizeof(struct unit_system));
@@ -6977,6 +6982,10 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
e->output_list_stf = output_list_stf;
}
+#ifdef EOS_PLANETARY
+ eos_init(&eos, e->physical_constants, e->snapshot_units, e->parameter_file);
+#endif
+
/* Want to force a rebuild before using this engine. Wait to repartition.*/
e->forcerebuild = 1;
e->forcerepart = 0;
diff --git a/src/engine.h b/src/engine.h
index aeb57c65ac36ff5ddbf4b74185adeb94f3d460da..cfd656712bdea84484101cf7e83795f795200f5f 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -223,6 +223,7 @@ struct engine {
char snapshot_base_name[PARSER_MAX_LINE_SIZE];
int snapshot_compression;
+ int snapshot_label_first;
int snapshot_label_delta;
struct unit_system *snapshot_units;
int snapshot_output_count;
diff --git a/src/equation_of_state/planetary/aneos.h b/src/equation_of_state/planetary/aneos.h
deleted file mode 100644
index 904288b2fdf3ba825cdc7d114ebb61cd42de198d..0000000000000000000000000000000000000000
--- a/src/equation_of_state/planetary/aneos.h
+++ /dev/null
@@ -1,144 +0,0 @@
-/*******************************************************************************
- * 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
index 61e23dc0b4eb82e9ae5c0869f7a10dfff97fc45e..644167bb4795a2a3d0fefe130cba93c64f29941b 100644
--- a/src/equation_of_state/planetary/equation_of_state.h
+++ b/src/equation_of_state/planetary/equation_of_state.h
@@ -39,6 +39,7 @@
#include "common_io.h"
#include "inline.h"
#include "physical_constants.h"
+#include "restart.h"
#include "units.h"
extern struct eos_parameters eos;
@@ -50,10 +51,15 @@ extern struct eos_parameters eos;
* @brief Master type for the planetary equation of state.
*/
enum eos_planetary_type_id {
+
+ /*! Tillotson */
eos_planetary_type_Til = 1,
+
+ /*! Hubbard & MacFarlane (1980) Uranus/Neptune */
eos_planetary_type_HM80 = 2,
- eos_planetary_type_ANEOS = 3,
- eos_planetary_type_SESAME = 4,
+
+ /*! SESAME */
+ eos_planetary_type_SESAME = 3,
};
/**
@@ -89,25 +95,26 @@ enum eos_planetary_material_id {
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,
-
- /*! MANEOS forsterite */
- eos_planetary_id_MANEOS_forsterite =
- eos_planetary_type_ANEOS * eos_planetary_type_factor + 1,
-
/* SESAME */
- /*! SESAME iron */
+ /*! SESAME iron 2140 */
eos_planetary_id_SESAME_iron =
eos_planetary_type_SESAME * eos_planetary_type_factor,
+
+ /*! SESAME basalt 7530 */
+ eos_planetary_id_SESAME_basalt =
+ eos_planetary_type_SESAME * eos_planetary_type_factor + 1,
+
+ /*! SESAME water 7154 */
+ eos_planetary_id_SESAME_water =
+ eos_planetary_type_SESAME * eos_planetary_type_factor + 2,
+
+ /*! Senft & Stewart (2008) SESAME-like water */
+ eos_planetary_id_SS08_water =
+ eos_planetary_type_SESAME * eos_planetary_type_factor + 3,
};
/* Individual EOS function headers. */
-#include "aneos.h"
#include "hm80.h"
#include "sesame.h"
#include "tillotson.h"
@@ -118,8 +125,7 @@ enum eos_planetary_material_id {
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;
+ struct SESAME_params SESAME_iron, SESAME_basalt, SESAME_water, SS08_water;
};
/**
@@ -190,35 +196,29 @@ gas_internal_energy_from_entropy(float density, float entropy,
};
break;
- /* ANEOS EoS */
- case eos_planetary_type_ANEOS:
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
/* Select the material */
switch (mat_id) {
- case eos_planetary_id_ANEOS_iron:
- return ANEOS_internal_energy_from_entropy(density, entropy,
- &eos.ANEOS_iron);
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_internal_energy_from_entropy(density, entropy,
+ &eos.SESAME_iron);
break;
- case eos_planetary_id_MANEOS_forsterite:
- return ANEOS_internal_energy_from_entropy(density, entropy,
- &eos.MANEOS_forsterite);
+ case eos_planetary_id_SESAME_basalt:
+ return SESAME_internal_energy_from_entropy(density, entropy,
+ &eos.SESAME_basalt);
break;
- default:
- error("Unknown material ID! mat_id = %d", mat_id);
- return 0.f;
- };
- break;
-
- /* SESAME EoS */
- case eos_planetary_type_SESAME:;
+ case eos_planetary_id_SESAME_water:
+ return SESAME_internal_energy_from_entropy(density, entropy,
+ &eos.SESAME_water);
+ break;
- /* Select the material */
- switch (mat_id) {
- case eos_planetary_id_SESAME_iron:
+ case eos_planetary_id_SS08_water:
return SESAME_internal_energy_from_entropy(density, entropy,
- &eos.SESAME_iron);
+ &eos.SS08_water);
break;
default:
@@ -294,34 +294,29 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
};
break;
- /* ANEOS EoS */
- case eos_planetary_type_ANEOS:
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
/* Select the material */
switch (mat_id) {
- case eos_planetary_id_ANEOS_iron:
- return ANEOS_pressure_from_entropy(density, entropy, &eos.ANEOS_iron);
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_pressure_from_entropy(density, entropy,
+ &eos.SESAME_iron);
break;
- case eos_planetary_id_MANEOS_forsterite:
- return ANEOS_pressure_from_entropy(density, entropy,
- &eos.MANEOS_forsterite);
+ case eos_planetary_id_SESAME_basalt:
+ return SESAME_pressure_from_entropy(density, entropy,
+ &eos.SESAME_basalt);
break;
- default:
- error("Unknown material ID! mat_id = %d", mat_id);
- return 0.f;
- };
- break;
-
- /* SESAME EoS */
- case eos_planetary_type_SESAME:;
+ case eos_planetary_id_SESAME_water:
+ return SESAME_pressure_from_entropy(density, entropy,
+ &eos.SESAME_water);
- /* Select the material */
- switch (mat_id) {
- case eos_planetary_id_SESAME_iron:
+ case eos_planetary_id_SS08_water:
return SESAME_pressure_from_entropy(density, entropy,
- &eos.SESAME_iron);
+ &eos.SS08_water);
+ break;
break;
default:
@@ -398,33 +393,25 @@ __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
};
break;
- /* ANEOS EoS */
- case eos_planetary_type_ANEOS:
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
/* Select the material */
switch (mat_id) {
- case eos_planetary_id_ANEOS_iron:
- return ANEOS_entropy_from_pressure(density, P, &eos.ANEOS_iron);
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_entropy_from_pressure(density, P, &eos.SESAME_iron);
break;
- case eos_planetary_id_MANEOS_forsterite:
- return ANEOS_entropy_from_pressure(density, P,
- &eos.MANEOS_forsterite);
+ case eos_planetary_id_SESAME_basalt:
+ return SESAME_entropy_from_pressure(density, P, &eos.SESAME_basalt);
break;
- default:
- error("Unknown material ID! mat_id = %d", mat_id);
- return 0.f;
- };
- break;
-
- /* SESAME EoS */
- case eos_planetary_type_SESAME:;
+ case eos_planetary_id_SESAME_water:
+ return SESAME_entropy_from_pressure(density, P, &eos.SESAME_water);
+ break;
- /* Select the material */
- switch (mat_id) {
- case eos_planetary_id_SESAME_iron:
- return SESAME_entropy_from_pressure(density, P, &eos.SESAME_iron);
+ case eos_planetary_id_SS08_water:
+ return SESAME_entropy_from_pressure(density, P, &eos.SS08_water);
break;
default:
@@ -501,35 +488,29 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
};
break;
- /* ANEOS EoS */
- case eos_planetary_type_ANEOS:
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
/* Select the material */
switch (mat_id) {
- case eos_planetary_id_ANEOS_iron:
- return ANEOS_soundspeed_from_entropy(density, entropy,
- &eos.ANEOS_iron);
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_soundspeed_from_entropy(density, entropy,
+ &eos.SESAME_iron);
break;
- case eos_planetary_id_MANEOS_forsterite:
- return ANEOS_soundspeed_from_entropy(density, entropy,
- &eos.MANEOS_forsterite);
+ case eos_planetary_id_SESAME_basalt:
+ return SESAME_soundspeed_from_entropy(density, entropy,
+ &eos.SESAME_basalt);
break;
- default:
- error("Unknown material ID! mat_id = %d", mat_id);
- return 0.f;
- };
- break;
-
- /* SESAME EoS */
- case eos_planetary_type_SESAME:;
+ case eos_planetary_id_SESAME_water:
+ return SESAME_soundspeed_from_entropy(density, entropy,
+ &eos.SESAME_water);
+ break;
- /* Select the material */
- switch (mat_id) {
- case eos_planetary_id_SESAME_iron:
+ case eos_planetary_id_SS08_water:
return SESAME_soundspeed_from_entropy(density, entropy,
- &eos.SESAME_iron);
+ &eos.SS08_water);
break;
default:
@@ -605,35 +586,29 @@ gas_entropy_from_internal_energy(float density, float u,
};
break;
- /* ANEOS EoS */
- case eos_planetary_type_ANEOS:
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
/* Select the material */
switch (mat_id) {
- case eos_planetary_id_ANEOS_iron:
- return ANEOS_entropy_from_internal_energy(density, u,
- &eos.ANEOS_iron);
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_entropy_from_internal_energy(density, u,
+ &eos.SESAME_iron);
break;
- case eos_planetary_id_MANEOS_forsterite:
- return ANEOS_entropy_from_internal_energy(density, u,
- &eos.MANEOS_forsterite);
+ case eos_planetary_id_SESAME_basalt:
+ return SESAME_entropy_from_internal_energy(density, u,
+ &eos.SESAME_basalt);
break;
- default:
- error("Unknown material ID! mat_id = %d", mat_id);
- return 0.f;
- };
- break;
-
- /* SESAME EoS */
- case eos_planetary_type_SESAME:;
+ case eos_planetary_id_SESAME_water:
+ return SESAME_entropy_from_internal_energy(density, u,
+ &eos.SESAME_water);
+ break;
- /* Select the material */
- switch (mat_id) {
- case eos_planetary_id_SESAME_iron:
+ case eos_planetary_id_SS08_water:
return SESAME_entropy_from_internal_energy(density, u,
- &eos.SESAME_iron);
+ &eos.SS08_water);
break;
default:
@@ -711,35 +686,29 @@ gas_pressure_from_internal_energy(float density, float u,
};
break;
- /* ANEOS EoS */
- case eos_planetary_type_ANEOS:
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
/* Select the material */
switch (mat_id) {
- case eos_planetary_id_ANEOS_iron:
- return ANEOS_pressure_from_internal_energy(density, u,
- &eos.ANEOS_iron);
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_pressure_from_internal_energy(density, u,
+ &eos.SESAME_iron);
break;
- case eos_planetary_id_MANEOS_forsterite:
- return ANEOS_pressure_from_internal_energy(density, u,
- &eos.MANEOS_forsterite);
+ case eos_planetary_id_SESAME_basalt:
+ return SESAME_pressure_from_internal_energy(density, u,
+ &eos.SESAME_basalt);
break;
- default:
- error("Unknown material ID! mat_id = %d", mat_id);
- return 0.f;
- };
- break;
-
- /* SESAME EoS */
- case eos_planetary_type_SESAME:;
+ case eos_planetary_id_SESAME_water:
+ return SESAME_pressure_from_internal_energy(density, u,
+ &eos.SESAME_water);
+ break;
- /* Select the material */
- switch (mat_id) {
- case eos_planetary_id_SESAME_iron:
+ case eos_planetary_id_SS08_water:
return SESAME_pressure_from_internal_energy(density, u,
- &eos.SESAME_iron);
+ &eos.SS08_water);
break;
default:
@@ -820,35 +789,29 @@ gas_internal_energy_from_pressure(float density, float P,
};
break;
- /* ANEOS EoS */
- case eos_planetary_type_ANEOS:
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
/* Select the material */
switch (mat_id) {
- case eos_planetary_id_ANEOS_iron:
- return ANEOS_internal_energy_from_pressure(density, P,
- &eos.ANEOS_iron);
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_internal_energy_from_pressure(density, P,
+ &eos.SESAME_iron);
break;
- case eos_planetary_id_MANEOS_forsterite:
- return ANEOS_internal_energy_from_pressure(density, P,
- &eos.MANEOS_forsterite);
+ case eos_planetary_id_SESAME_basalt:
+ return SESAME_internal_energy_from_pressure(density, P,
+ &eos.SESAME_basalt);
break;
- default:
- error("Unknown material ID! mat_id = %d", mat_id);
- return 0.f;
- };
- break;
-
- /* SESAME EoS */
- case eos_planetary_type_SESAME:;
+ case eos_planetary_id_SESAME_water:
+ return SESAME_internal_energy_from_pressure(density, P,
+ &eos.SESAME_water);
+ break;
- /* Select the material */
- switch (mat_id) {
- case eos_planetary_id_SESAME_iron:
+ case eos_planetary_id_SS08_water:
return SESAME_internal_energy_from_pressure(density, P,
- &eos.SESAME_iron);
+ &eos.SS08_water);
break;
default:
@@ -930,35 +893,29 @@ gas_soundspeed_from_internal_energy(float density, float u,
};
break;
- /* ANEOS EoS */
- case eos_planetary_type_ANEOS:
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
/* Select the material */
switch (mat_id) {
- case eos_planetary_id_ANEOS_iron:
- return ANEOS_soundspeed_from_internal_energy(density, u,
- &eos.ANEOS_iron);
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_soundspeed_from_internal_energy(density, u,
+ &eos.SESAME_iron);
break;
- case eos_planetary_id_MANEOS_forsterite:
- return ANEOS_soundspeed_from_internal_energy(density, u,
- &eos.MANEOS_forsterite);
+ case eos_planetary_id_SESAME_basalt:
+ return SESAME_soundspeed_from_internal_energy(density, u,
+ &eos.SESAME_basalt);
break;
- default:
- error("Unknown material ID! mat_id = %d", mat_id);
- return 0.f;
- };
- break;
-
- /* SESAME EoS */
- case eos_planetary_type_SESAME:;
+ case eos_planetary_id_SESAME_water:
+ return SESAME_soundspeed_from_internal_energy(density, u,
+ &eos.SESAME_water);
+ break;
- /* Select the material */
- switch (mat_id) {
- case eos_planetary_id_SESAME_iron:
+ case eos_planetary_id_SS08_water:
return SESAME_soundspeed_from_internal_energy(density, u,
- &eos.SESAME_iron);
+ &eos.SS08_water);
break;
default:
@@ -1034,33 +991,26 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
};
break;
- /* ANEOS EoS */
- case eos_planetary_type_ANEOS:
+ /* SESAME EoS */
+ case eos_planetary_type_SESAME:;
/* Select the material */
switch (mat_id) {
- case eos_planetary_id_ANEOS_iron:
- return ANEOS_soundspeed_from_pressure(density, P, &eos.ANEOS_iron);
+ case eos_planetary_id_SESAME_iron:
+ return SESAME_soundspeed_from_pressure(density, P, &eos.SESAME_iron);
break;
- case eos_planetary_id_MANEOS_forsterite:
- return ANEOS_soundspeed_from_pressure(density, P,
- &eos.MANEOS_forsterite);
+ case eos_planetary_id_SESAME_basalt:
+ return SESAME_soundspeed_from_pressure(density, P,
+ &eos.SESAME_basalt);
break;
- default:
- error("Unknown material ID! mat_id = %d", mat_id);
- return 0.f;
- };
- break;
-
- /* SESAME EoS */
- case eos_planetary_type_SESAME:;
+ case eos_planetary_id_SESAME_water:
+ return SESAME_soundspeed_from_pressure(density, P, &eos.SESAME_water);
+ break;
- /* Select the material */
- switch (mat_id) {
- case eos_planetary_id_SESAME_iron:
- return SESAME_soundspeed_from_pressure(density, P, &eos.SESAME_iron);
+ case eos_planetary_id_SS08_water:
+ return SESAME_soundspeed_from_pressure(density, P, &eos.SS08_water);
break;
default:
@@ -1089,6 +1039,10 @@ __attribute__((always_inline)) INLINE static void eos_init(
char HM80_HHe_table_file[PARSER_MAX_LINE_SIZE];
char HM80_ice_table_file[PARSER_MAX_LINE_SIZE];
char HM80_rock_table_file[PARSER_MAX_LINE_SIZE];
+ char SESAME_iron_table_file[PARSER_MAX_LINE_SIZE];
+ char SESAME_basalt_table_file[PARSER_MAX_LINE_SIZE];
+ char SESAME_water_table_file[PARSER_MAX_LINE_SIZE];
+ char SS08_water_table_file[PARSER_MAX_LINE_SIZE];
// Set the parameters and material IDs, load tables, etc. for each material
// and convert to internal units
@@ -1116,30 +1070,49 @@ __attribute__((always_inline)) INLINE static void eos_init(
parser_get_param_string(params, "EoS:planetary_HM80_rock_table_file",
HM80_rock_table_file);
- 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);
+ load_table_HM80(&e->HM80_HHe, HM80_HHe_table_file);
+ load_table_HM80(&e->HM80_ice, HM80_ice_table_file);
+ load_table_HM80(&e->HM80_rock, HM80_rock_table_file);
+
+ prepare_table_HM80(&e->HM80_HHe);
+ prepare_table_HM80(&e->HM80_ice);
+ prepare_table_HM80(&e->HM80_rock);
convert_units_HM80(&e->HM80_HHe, us);
convert_units_HM80(&e->HM80_ice, us);
convert_units_HM80(&e->HM80_rock, us);
}
- // ANEOS
- if (parser_get_opt_param_int(params, "EoS:planetary_use_ANEOS", 0)) {
- set_ANEOS_iron(&e->ANEOS_iron, eos_planetary_id_ANEOS_iron);
- set_MANEOS_forsterite(&e->MANEOS_forsterite,
- eos_planetary_id_MANEOS_forsterite);
-
- convert_units_ANEOS(&e->ANEOS_iron, us);
- convert_units_ANEOS(&e->MANEOS_forsterite, us);
- }
-
// SESAME
if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME", 0)) {
set_SESAME_iron(&e->SESAME_iron, eos_planetary_id_SESAME_iron);
+ set_SESAME_basalt(&e->SESAME_basalt, eos_planetary_id_SESAME_basalt);
+ set_SESAME_water(&e->SESAME_water, eos_planetary_id_SESAME_water);
+ set_SS08_water(&e->SESAME_water, eos_planetary_id_SS08_water);
+
+ parser_get_param_string(params, "EoS:planetary_SESAME_iron_table_file",
+ SESAME_iron_table_file);
+ parser_get_param_string(params, "EoS:planetary_SESAME_basalt_table_file",
+ SESAME_basalt_table_file);
+ parser_get_param_string(params, "EoS:planetary_SESAME_water_table_file",
+ SESAME_water_table_file);
+ parser_get_param_string(params, "EoS:planetary_SS08_water_table_file",
+ SS08_water_table_file);
+
+ load_table_SESAME(&e->SESAME_iron, SESAME_iron_table_file);
+ load_table_SESAME(&e->SESAME_basalt, SESAME_basalt_table_file);
+ load_table_SESAME(&e->SESAME_water, SESAME_water_table_file);
+ load_table_SESAME(&e->SS08_water, SS08_water_table_file);
+
+ prepare_table_SESAME(&e->SESAME_iron);
+ prepare_table_SESAME(&e->SESAME_basalt);
+ prepare_table_SESAME(&e->SESAME_water);
+ prepare_table_SESAME(&e->SS08_water);
convert_units_SESAME(&e->SESAME_iron, us);
+ convert_units_SESAME(&e->SESAME_basalt, us);
+ convert_units_SESAME(&e->SESAME_water, us);
+ convert_units_SESAME(&e->SS08_water, us);
}
}
diff --git a/src/equation_of_state/planetary/get_eos_tables.sh b/src/equation_of_state/planetary/get_eos_tables.sh
new file mode 100755
index 0000000000000000000000000000000000000000..c0a751252bb060341a01ac70320a16251069a84e
--- /dev/null
+++ b/src/equation_of_state/planetary/get_eos_tables.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# Download the tables of the publicly available planetary equations of state
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/planetary_HM80_HHe.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/planetary_HM80_ice.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/planetary_HM80_rock.txt
+
+mv planetary_HM80_HHe.txt ../../../examples/
+mv planetary_HM80_ice.txt ../../../examples/
+mv planetary_HM80_rock.txt ../../../examples/
diff --git a/src/equation_of_state/planetary/hm80.h b/src/equation_of_state/planetary/hm80.h
index 0131bab6c447e5a8898e29e13dc3f8f6e1c897c6..3f9600709836b99c74a4857653214497e0aaa773 100644
--- a/src/equation_of_state/planetary/hm80.h
+++ b/src/equation_of_state/planetary/hm80.h
@@ -41,75 +41,79 @@
// Hubbard & MacFarlane (1980) parameters
struct HM80_params {
- float *table_P_rho_u;
+ float *table_log_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;
+ log_u_max, log_u_step, inv_log_u_step, bulk_mod, P_min_for_c_min;
enum eos_planetary_material_id mat_id;
};
-// Parameter values for each material (cgs units)
+// Parameter values for each material (SI 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;
+ mat->bulk_mod = 0.f;
+ mat->P_min_for_c_min = 1e3f;
}
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;
+ mat->bulk_mod = 2.0e9f;
+ mat->P_min_for_c_min = 0.f;
}
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;
+ mat->bulk_mod = 3.49e10f;
+ mat->P_min_for_c_min = 0.f;
}
// 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 * mat->num_u * sizeof(float *));
+INLINE static void load_table_HM80(struct HM80_params *mat, char *table_file) {
+
+ /* File contents:
+ header (four lines)
+ log_rho_min log_rho_max num_rho log_u_min log_u_max num_u (SI)
+ P_0_0 P_0_1 ... P_0_num_u # Array of pressures (Pa)
+ P_1_0 ... ... P_1_num_u
+ ... ... ... ...
+ P_num_rho_0 ... P_num_rho_num_u
+ T_0_0 T_0_1 ... T_0_num_u # Array of temperatures (K)
+ T_1_0 ... ... T_1_num_u
+ ... ... ... ...
+ T_num_rho_0 ... T_num_rho_num_u
+ */
// 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 * mat->num_rho + j]);
+
+ // Ignore header lines
+ char buffer[100];
+ for (int i = 0; i < 4; i++) {
+ fgets(buffer, 100, f);
+ }
+
+ // Table properties
+ c = fscanf(f, "%f %f %d %f %f %d", &mat->log_rho_min, &mat->log_rho_max,
+ &mat->num_rho, &mat->log_u_min, &mat->log_u_max, &mat->num_u);
+ if (c != 6) {
+ error("Failed to read EOS table %s", table_file);
+ }
+ mat->log_rho_step =
+ (mat->log_rho_max - mat->log_rho_min) / (mat->num_rho - 1);
+ mat->log_u_step = (mat->log_u_max - mat->log_u_min) / (mat->num_u - 1);
+ mat->inv_log_rho_step = 1.f / mat->log_rho_step;
+ mat->inv_log_u_step = 1.f / mat->log_u_step;
+
+ // Allocate table memory
+ mat->table_log_P_rho_u =
+ (float *)malloc(mat->num_rho * mat->num_u * sizeof(float));
+
+ // Pressures (not log yet)
+ for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
+ for (int i_u = 0; i_u < mat->num_u; i_u++) {
+ c = fscanf(f, "%f", &mat->table_log_P_rho_u[i_rho * mat->num_u + i_u]);
if (c != 1) {
error("Failed to read EOS table");
}
@@ -118,33 +122,49 @@ INLINE static void load_HM80_table(struct HM80_params *mat, char *table_file) {
fclose(f);
}
-// Convert from cgs to internal units
+// Misc. modifications
+INLINE static void prepare_table_HM80(struct HM80_params *mat) {
+
+ // Convert pressures to log(pressure)
+ for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
+ for (int i_u = 0; i_u < mat->num_u; i_u++) {
+ mat->table_log_P_rho_u[i_rho * mat->num_u + i_u] =
+ logf(mat->table_log_P_rho_u[i_rho * mat->num_u + i_u]);
+ }
+ }
+}
+
+// Convert 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
+ struct unit_system si;
+ units_init_si(&si);
- // 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));
+ // All table values in SI
+ mat->log_rho_min += logf(units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY) /
+ units_cgs_conversion_factor(us, UNIT_CONV_DENSITY));
+ mat->log_rho_max += logf(units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY) /
+ units_cgs_conversion_factor(us, UNIT_CONV_DENSITY));
- // Table energies in SI
mat->log_u_min +=
- logf(J_kg_to_erg_g /
+ logf(units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS) /
units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
mat->log_u_max +=
- logf(J_kg_to_erg_g /
+ logf(units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS) /
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 * mat->num_rho + j] *=
- Mbar_to_Ba / units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
+ for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
+ for (int i_u = 0; i_u < mat->num_u; i_u++) {
+ mat->table_log_P_rho_u[i_rho * mat->num_u + i_u] +=
+ logf(units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE) /
+ units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE));
}
}
- mat->bulk_mod /= units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
+ mat->bulk_mod *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE) /
+ units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
+ mat->P_min_for_c_min *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE) /
+ units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
}
// gas_internal_energy_from_entropy
@@ -153,7 +173,7 @@ INLINE static float HM80_internal_energy_from_entropy(
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// gas_pressure_from_entropy
@@ -162,7 +182,7 @@ INLINE static float HM80_pressure_from_entropy(float density, float entropy,
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// gas_entropy_from_pressure
@@ -171,7 +191,7 @@ INLINE static float HM80_entropy_from_pressure(float density, float pressure,
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// gas_soundspeed_from_entropy
@@ -180,75 +200,62 @@ INLINE static float HM80_soundspeed_from_entropy(
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// gas_entropy_from_internal_energy
INLINE static float HM80_entropy_from_internal_energy(
float density, float u, const struct HM80_params *mat) {
- return 0;
+ return 0.f;
}
// gas_pressure_from_internal_energy
INLINE static float HM80_pressure_from_internal_energy(
float density, float u, const struct HM80_params *mat) {
- float P;
+ float log_P, log_P_1, log_P_2, log_P_3, log_P_4;
if (u <= 0.f) {
return 0.f;
}
- int rho_idx, u_idx;
+ int idx_rho, idx_u;
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);
+ // 2D interpolation (bilinear with log(rho), log(u)) to find P(rho, u)
+ idx_rho = floor((log_rho - mat->log_rho_min) * mat->inv_log_rho_step);
+ idx_u = floor((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[u_idx] +
- intp_u * mat->table_P_rho_u[u_idx + 1]) +
- log_rho - mat->log_rho_min);
- if (u_idx < 0) { // and too-low u
- P = 0.f;
- }
- } else if (u_idx < 0) { // Too-low u
- P = 0.f;
+ // If outside the table then extrapolate from the edge and edge-but-one values
+ if (idx_rho <= -1) {
+ idx_rho = 0;
+ } else if (idx_rho >= mat->num_rho - 1) {
+ idx_rho = mat->num_rho - 2;
}
- // 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 + mat->num_u - 1];
- } else {
- P = mat->table_P_rho_u[(mat->num_rho - 1) * mat->num_u + u_idx];
- }
- } else if (u_idx >= mat->num_u - 1) { // Too-high u
- P = mat->table_P_rho_u[rho_idx * mat->num_u + 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 * mat->num_u + u_idx] +
- intp_u * mat->table_P_rho_u[rho_idx * mat->num_u + u_idx + 1]) +
- intp_rho *
- ((1 - intp_u) *
- mat->table_P_rho_u[(rho_idx + 1) * mat->num_u + u_idx] +
- intp_u *
- mat->table_P_rho_u[(rho_idx + 1) * mat->num_u + u_idx + 1]);
+ if (idx_u <= -1) {
+ idx_u = 0;
+ } else if (idx_u >= mat->num_u - 1) {
+ idx_u = mat->num_u - 2;
}
- return P;
+ intp_rho = (log_rho - mat->log_rho_min - idx_rho * mat->log_rho_step) *
+ mat->inv_log_rho_step;
+ intp_u =
+ (log_u - mat->log_u_min - idx_u * mat->log_u_step) * mat->inv_log_u_step;
+
+ // Table values
+ log_P_1 = mat->table_log_P_rho_u[idx_rho * mat->num_u + idx_u];
+ log_P_2 = mat->table_log_P_rho_u[idx_rho * mat->num_u + idx_u + 1];
+ log_P_3 = mat->table_log_P_rho_u[(idx_rho + 1) * mat->num_u + idx_u];
+ log_P_4 = mat->table_log_P_rho_u[(idx_rho + 1) * mat->num_u + idx_u + 1];
+
+ log_P = (1.f - intp_rho) * ((1.f - intp_u) * log_P_1 + intp_u * log_P_2) +
+ intp_rho * ((1.f - intp_u) * log_P_3 + intp_u * log_P_4);
+
+ return expf(log_P);
}
// gas_internal_energy_from_pressure
@@ -257,7 +264,7 @@ INLINE static float HM80_internal_energy_from_pressure(
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// gas_soundspeed_from_internal_energy
@@ -274,6 +281,10 @@ INLINE static float HM80_soundspeed_from_internal_energy(
else {
P = HM80_pressure_from_internal_energy(density, u, mat);
c = sqrtf(hydro_gamma * P / density);
+
+ if (c <= 0) {
+ c = sqrtf(hydro_gamma * mat->P_min_for_c_min / density);
+ }
}
return c;
@@ -283,18 +294,9 @@ INLINE static float HM80_soundspeed_from_internal_energy(
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);
- }
+ error("This EOS function is not yet implemented!");
- return c;
+ return 0.f;
}
#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
index 76574c2ad00282a82649705cd8a2b5a1b428d867..4deaded44314d6b6204c25a18a6eacc4bdc6ac9a 100644
--- a/src/equation_of_state/planetary/sesame.h
+++ b/src/equation_of_state/planetary/sesame.h
@@ -40,21 +40,228 @@
#include "inline.h"
#include "physical_constants.h"
#include "units.h"
+#include "utilities.h"
// SESAME parameters
struct SESAME_params {
+ float *table_log_rho;
+ float *table_log_u_rho_T;
+ float *table_P_rho_T;
+ float *table_c_rho_T;
+ float *table_s_rho_T;
+ int num_rho, num_T;
+ float P_tiny, c_tiny;
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) {
+ // SESAME 2140
mat->mat_id = mat_id;
}
+INLINE static void set_SESAME_basalt(struct SESAME_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ // SESAME 7530
+ mat->mat_id = mat_id;
+}
+INLINE static void set_SESAME_water(struct SESAME_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ // SESAME 7154
+ mat->mat_id = mat_id;
+}
+INLINE static void set_SS08_water(struct SESAME_params *mat,
+ enum eos_planetary_material_id mat_id) {
+ // Senft & Stewart (2008)
+ mat->mat_id = mat_id;
+}
+
+// Read the tables from file
+INLINE static void load_table_SESAME(struct SESAME_params *mat,
+ char *table_file) {
+
+ // Load table contents from file
+ FILE *f = fopen(table_file, "r");
+ int c;
+
+ // Ignore header lines
+ char buffer[100];
+ for (int i = 0; i < 5; i++) {
+ fgets(buffer, 100, f);
+ }
+ float ignore;
+
+ // Table properties
+ c = fscanf(f, "%d %d", &mat->num_rho, &mat->num_T);
+ if (c != 2) {
+ error("Failed to read EOS table %s", table_file);
+ }
+
+ // Ignore the first elements of rho = 0, T = 0
+ mat->num_rho--;
+ mat->num_T--;
+
+ // Allocate table memory
+ mat->table_log_rho = (float *)malloc(mat->num_rho * sizeof(float));
+ mat->table_log_u_rho_T =
+ (float *)malloc(mat->num_rho * mat->num_T * sizeof(float));
+ mat->table_P_rho_T =
+ (float *)malloc(mat->num_rho * mat->num_T * sizeof(float));
+ mat->table_c_rho_T =
+ (float *)malloc(mat->num_rho * mat->num_T * sizeof(float));
+ mat->table_s_rho_T =
+ (float *)malloc(mat->num_rho * mat->num_T * sizeof(float));
+
+ // Densities (not log yet)
+ for (int i_rho = -1; i_rho < mat->num_rho; i_rho++) {
+ // Ignore the first elements of rho = 0, T = 0
+ if (i_rho == -1) {
+ c = fscanf(f, "%f", &ignore);
+ if (c != 1) {
+ error("Failed to read EOS table %s", table_file);
+ }
+ } else {
+ c = fscanf(f, "%f", &mat->table_log_rho[i_rho]);
+ if (c != 1) {
+ error("Failed to read EOS table %s", table_file);
+ }
+ }
+ }
+
+ // Temperatures (ignored)
+ for (int i_T = -1; i_T < mat->num_T; i_T++) {
+ c = fscanf(f, "%f", &ignore);
+ if (c != 1) {
+ error("Failed to read EOS table %s", table_file);
+ }
+ }
+
+ // Sp. int. energies (not log yet), pressures, sound speeds, and entropies
+ for (int i_T = -1; i_T < mat->num_T; i_T++) {
+ for (int i_rho = -1; i_rho < mat->num_rho; i_rho++) {
+ // Ignore the first elements of rho = 0, T = 0
+ if ((i_T == -1) || (i_rho == -1)) {
+ c = fscanf(f, "%f %f %f %f", &ignore, &ignore, &ignore, &ignore);
+ if (c != 4) {
+ error("Failed to read EOS table %s", table_file);
+ }
+ } else {
+ c = fscanf(f, "%f %f %f %f",
+ &mat->table_log_u_rho_T[i_rho * mat->num_T + i_T],
+ &mat->table_P_rho_T[i_rho * mat->num_T + i_T],
+ &mat->table_c_rho_T[i_rho * mat->num_T + i_T],
+ &mat->table_s_rho_T[i_rho * mat->num_T + i_T]);
+ if (c != 4) {
+ error("Failed to read EOS table %s", table_file);
+ }
+ }
+ }
+ }
+
+ fclose(f);
+}
-// Convert from cgs to internal units
+// Misc. modifications
+INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
+
+ // Convert densities to log(density)
+ for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
+ mat->table_log_rho[i_rho] = logf(mat->table_log_rho[i_rho]);
+ }
+
+ // Convert sp. int. energies to log(sp. int. energy)
+ for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
+ for (int i_T = 0; i_T < mat->num_T; i_T++) {
+ // If not positive then set very small for the log
+ if (mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] <= 0) {
+ mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] = 1.f;
+ }
+
+ mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] =
+ logf(mat->table_log_u_rho_T[i_rho * mat->num_T + i_T]);
+ }
+ }
+
+ // Tiny pressure and soundspeed, initialise in the middle
+ mat->P_tiny =
+ mat->table_P_rho_T[mat->num_rho / 2 * mat->num_T + mat->num_T / 2];
+ mat->c_tiny =
+ mat->table_c_rho_T[mat->num_rho / 2 * mat->num_T + mat->num_T / 2];
+
+ // Enforce that the 1D arrays of u (at each rho) are monotonic
+ // This is necessary because, for some high-density u slices at very low T,
+ // u decreases (very slightly) with T, which makes the interpolation fail
+ for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
+ for (int i_T = mat->num_T - 1; i_T > 0; i_T--) {
+
+ // If the one-lower-T u is greater than this u
+ if (mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] <
+ mat->table_log_u_rho_T[i_rho * mat->num_T + i_T - 1]) {
+
+ // Replace it and all elements below it with that value
+ for (int j_u = 0; j_u < i_T; j_u++) {
+ mat->table_log_u_rho_T[i_rho * mat->num_T + j_u] =
+ mat->table_log_u_rho_T[i_rho * mat->num_T + i_T];
+ }
+ break;
+ }
+
+ // Smallest positive pressure and sound speed
+ if ((mat->table_P_rho_T[i_rho * mat->num_T + i_T] < mat->P_tiny) &&
+ (mat->table_P_rho_T[i_rho * mat->num_T + i_T] > 0)) {
+ mat->P_tiny = mat->table_P_rho_T[i_rho * mat->num_T + i_T];
+ }
+ if ((mat->table_c_rho_T[i_rho * mat->num_T + i_T] < mat->c_tiny) &&
+ (mat->table_c_rho_T[i_rho * mat->num_T + i_T] > 0)) {
+ mat->c_tiny = mat->table_c_rho_T[i_rho * mat->num_T + i_T];
+ }
+ }
+ }
+
+ // Tiny pressure to allow interpolation near non-positive values
+ mat->P_tiny *= 1e-3f;
+ mat->c_tiny *= 1e-3f;
+}
+
+// Convert to internal units
INLINE static void convert_units_SESAME(struct SESAME_params *mat,
- const struct unit_system *us) {}
+ const struct unit_system *us) {
+
+ struct unit_system si;
+ units_init_si(&si);
+
+ // All table values in SI
+ // Densities (log)
+ for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
+ mat->table_log_rho[i_rho] +=
+ logf(units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY) /
+ units_cgs_conversion_factor(us, UNIT_CONV_DENSITY));
+ }
+
+ // Sp. Int. Energies (log), pressures, and sound speeds
+ for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
+ for (int i_T = 0; i_T < mat->num_T; i_T++) {
+ mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] += logf(
+ units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS) /
+ units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
+ mat->table_P_rho_T[i_rho * mat->num_T + i_T] *=
+ units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE) /
+ units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
+ mat->table_c_rho_T[i_rho * mat->num_T + i_T] *=
+ units_cgs_conversion_factor(&si, UNIT_CONV_SPEED) /
+ units_cgs_conversion_factor(us, UNIT_CONV_SPEED);
+ mat->table_s_rho_T[i_rho * mat->num_T + i_T] *=
+ units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS) /
+ units_cgs_conversion_factor(us, UNIT_CONV_ENTROPY);
+ }
+ }
+
+ // Tiny pressure and sound speed
+ mat->P_tiny *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE) /
+ units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
+ mat->c_tiny *= units_cgs_conversion_factor(&si, UNIT_CONV_SPEED) /
+ units_cgs_conversion_factor(us, UNIT_CONV_SPEED);
+}
// gas_internal_energy_from_entropy
INLINE static float SESAME_internal_energy_from_entropy(
@@ -62,7 +269,7 @@ INLINE static float SESAME_internal_energy_from_entropy(
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// gas_pressure_from_entropy
@@ -71,7 +278,7 @@ INLINE static float SESAME_pressure_from_entropy(
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// gas_entropy_from_pressure
@@ -80,7 +287,7 @@ INLINE static float SESAME_entropy_from_pressure(
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// gas_soundspeed_from_entropy
@@ -89,25 +296,109 @@ INLINE static float SESAME_soundspeed_from_entropy(
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// 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;
+ return 0.f;
}
// 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;
+ float P, P_1, P_2, P_3, P_4;
+
+ if (u <= 0.f) {
+ return 0.f;
+ }
+
+ int idx_rho, idx_u_1, idx_u_2;
+ float intp_rho, intp_u_1, intp_u_2;
+ const float log_rho = logf(density);
+ const float log_u = logf(u);
+
+ // 2D interpolation (bilinear with log(rho), log(u)) to find P(rho, u)
+ // Density index
+ idx_rho =
+ find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
+
+ // Sp. int. energy at this and the next density (in relevant slice of u array)
+ idx_u_1 = find_value_in_monot_incr_array(
+ log_u, mat->table_log_u_rho_T + idx_rho * mat->num_T, mat->num_T);
+ idx_u_2 = find_value_in_monot_incr_array(
+ log_u, mat->table_log_u_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
+
+ // If outside the table then extrapolate from the edge and edge-but-one values
+ if (idx_rho <= -1) {
+ idx_rho = 0;
+ } else if (idx_rho >= mat->num_rho) {
+ idx_rho = mat->num_rho - 2;
+ }
+ if (idx_u_1 <= -1) {
+ idx_u_1 = 0;
+ } else if (idx_u_1 >= mat->num_T) {
+ idx_u_1 = mat->num_T - 2;
+ }
+ if (idx_u_2 <= -1) {
+ idx_u_2 = 0;
+ } else if (idx_u_2 >= mat->num_T) {
+ idx_u_2 = mat->num_T - 2;
+ }
+
+ intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
+ (mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]);
+ intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) /
+ (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] -
+ mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]);
+ intp_u_2 =
+ (log_u - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) /
+ (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] -
+ mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]);
+
+ // Table values
+ P_1 = mat->table_P_rho_T[idx_rho * mat->num_T + idx_u_1];
+ P_2 = mat->table_P_rho_T[idx_rho * mat->num_T + idx_u_1 + 1];
+ P_3 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
+ P_4 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
+
+ // If more than two table values are non-positive then return zero
+ int num_non_pos = 0;
+ if (P_1 <= 0.f) num_non_pos++;
+ if (P_2 <= 0.f) num_non_pos++;
+ if (P_3 <= 0.f) num_non_pos++;
+ if (P_4 <= 0.f) num_non_pos++;
+ if (num_non_pos > 2) {
+ return 0.f;
+ }
+ // If just one or two are non-positive then replace them with a tiny value
+ else if (num_non_pos > 0) {
+ // Unless already trying to extrapolate in which case return zero
+ if ((intp_rho < 0.f) || (intp_u_1 < 0.f) || (intp_u_2 < 0.f)) {
+ return 0.f;
+ }
+ if (P_1 <= 0.f) P_1 = mat->P_tiny;
+ if (P_2 <= 0.f) P_2 = mat->P_tiny;
+ if (P_3 <= 0.f) P_3 = mat->P_tiny;
+ if (P_4 <= 0.f) P_4 = mat->P_tiny;
+ }
+
+ // Interpolate with the log values
+ P_1 = logf(P_1);
+ P_2 = logf(P_2);
+ P_3 = logf(P_3);
+ P_4 = logf(P_4);
+
+ P = (1.f - intp_rho) * ((1.f - intp_u_1) * P_1 + intp_u_1 * P_2) +
+ intp_rho * ((1.f - intp_u_2) * P_3 + intp_u_2 * P_4);
+
+ // Convert back from log
+ P = expf(P);
+
+ return P;
}
// gas_internal_energy_from_pressure
@@ -116,16 +407,102 @@ INLINE static float SESAME_internal_energy_from_pressure(
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// 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;
+ float c, c_1, c_2, c_3, c_4;
+
+ if (u <= 0.f) {
+ return 0.f;
+ }
+
+ int idx_rho, idx_u_1, idx_u_2;
+ float intp_rho, intp_u_1, intp_u_2;
+ const float log_rho = logf(density);
+ const float log_u = logf(u);
+
+ // 2D interpolation (bilinear with log(rho), log(u)) to find c(rho, u)
+ // Density index
+ idx_rho =
+ find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
+
+ // Sp. int. energy at this and the next density (in relevant slice of u array)
+ idx_u_1 = find_value_in_monot_incr_array(
+ log_u, mat->table_log_u_rho_T + idx_rho * mat->num_T, mat->num_T);
+ idx_u_2 = find_value_in_monot_incr_array(
+ log_u, mat->table_log_u_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
+
+ // If outside the table then extrapolate from the edge and edge-but-one values
+ if (idx_rho <= -1) {
+ idx_rho = 0;
+ } else if (idx_rho >= mat->num_rho) {
+ idx_rho = mat->num_rho - 2;
+ }
+ if (idx_u_1 <= -1) {
+ idx_u_1 = 0;
+ } else if (idx_u_1 >= mat->num_T) {
+ idx_u_1 = mat->num_T - 2;
+ }
+ if (idx_u_2 <= -1) {
+ idx_u_2 = 0;
+ } else if (idx_u_2 >= mat->num_T) {
+ idx_u_2 = mat->num_T - 2;
+ }
+
+ intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
+ (mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]);
+ intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) /
+ (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] -
+ mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]);
+ intp_u_2 =
+ (log_u - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) /
+ (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] -
+ mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]);
+
+ // Table values
+ c_1 = mat->table_c_rho_T[idx_rho * mat->num_T + idx_u_1];
+ c_2 = mat->table_c_rho_T[idx_rho * mat->num_T + idx_u_1 + 1];
+ c_3 = mat->table_c_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
+ c_4 = mat->table_c_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
+
+ // If more than two table values are non-positive then return zero
+ int num_non_pos = 0;
+ if (c_1 <= 0.f) num_non_pos++;
+ if (c_2 <= 0.f) num_non_pos++;
+ if (c_3 <= 0.f) num_non_pos++;
+ if (c_4 <= 0.f) num_non_pos++;
+ if (num_non_pos > 2) {
+ return mat->c_tiny;
+ }
+ // If just one or two are non-positive then replace them with a tiny value
+ else if (num_non_pos > 0) {
+ // Unless already trying to extrapolate in which case return zero
+ if ((intp_rho < 0.f) || (intp_u_1 < 0.f) || (intp_u_2 < 0.f)) {
+ return mat->c_tiny;
+ }
+ if (c_1 <= 0.f) c_1 = mat->c_tiny;
+ if (c_2 <= 0.f) c_2 = mat->c_tiny;
+ if (c_3 <= 0.f) c_3 = mat->c_tiny;
+ if (c_4 <= 0.f) c_4 = mat->c_tiny;
+ }
+
+ // Interpolate with the log values
+ c_1 = logf(c_1);
+ c_2 = logf(c_2);
+ c_3 = logf(c_3);
+ c_4 = logf(c_4);
+
+ c = (1.f - intp_rho) * ((1.f - intp_u_1) * c_1 + intp_u_1 * c_2) +
+ intp_rho * ((1.f - intp_u_2) * c_3 + intp_u_2 * c_4);
+
+ // Convert back from log
+ c = expf(c);
+
+ return c;
}
// gas_soundspeed_from_pressure
@@ -134,7 +511,7 @@ INLINE static float SESAME_soundspeed_from_pressure(
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
#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
index d5b6d5c35d5edf9e114fe7f010c4f5b1e2327a83..1a4210699380b3b0398506dde7fce6ca8055e4dc 100644
--- a/src/equation_of_state/planetary/tillotson.h
+++ b/src/equation_of_state/planetary/tillotson.h
@@ -41,22 +41,22 @@
// Tillotson parameters
struct Til_params {
- float rho_0, a, b, A, B, E_0, E_iv, E_cv, alpha, beta, eta_min, P_min;
+ float rho_0, a, b, A, B, u_0, u_iv, u_cv, alpha, beta, eta_min, P_min;
enum eos_planetary_material_id mat_id;
};
-// Parameter values for each material (cgs units)
+// Parameter values for each material (SI 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->rho_0 = 7800.0f;
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->A = 1.28e11f;
+ mat->B = 1.05e11f;
+ mat->u_0 = 9.5e9f;
+ mat->u_iv = 2.4e9f;
+ mat->u_cv = 8.67e9f;
mat->alpha = 5.0f;
mat->beta = 5.0f;
mat->eta_min = 0.0f;
@@ -65,14 +65,14 @@ INLINE static void set_Til_iron(struct Til_params *mat,
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->rho_0 = 2680.0f;
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->A = 1.8e10f;
+ mat->B = 1.8e10f;
+ mat->u_0 = 1.6e10f;
+ mat->u_iv = 3.5e9f;
+ mat->u_cv = 1.8e10f;
mat->alpha = 5.0f;
mat->beta = 5.0f;
mat->eta_min = 0.0f;
@@ -81,30 +81,43 @@ INLINE static void set_Til_granite(struct Til_params *mat,
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->rho_0 = 998.0f;
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->A = 2.18e9f;
+ mat->B = 1.325e10f;
+ mat->u_0 = 7.0e9f;
+ mat->u_iv = 4.19e8f;
+ mat->u_cv = 2.69e9f;
mat->alpha = 10.0f;
mat->beta = 5.0f;
- mat->eta_min = 0.915f;
+ mat->eta_min = 0.9f;
mat->P_min = 0.0f;
}
-// Convert from cgs to internal units
+// Convert to internal units
INLINE static void convert_units_Til(struct Til_params *mat,
const struct unit_system *us) {
+ struct unit_system si;
+ units_init_si(&si);
+
+ // SI to cgs
+ mat->rho_0 *= units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY);
+ mat->A *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE);
+ mat->B *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE);
+ mat->u_0 *= units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ mat->u_iv *= units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ mat->u_cv *= units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ mat->P_min *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE);
+
+ // cgs to internal
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->u_0 /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ mat->u_iv /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ mat->u_cv /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
mat->P_min /= units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
}
@@ -114,7 +127,7 @@ INLINE static float Til_internal_energy_from_entropy(
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// gas_pressure_from_entropy
@@ -123,7 +136,7 @@ INLINE static float Til_pressure_from_entropy(float density, float entropy,
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// gas_entropy_from_pressure
@@ -132,7 +145,7 @@ INLINE static float Til_entropy_from_pressure(float density, float pressure,
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// gas_soundspeed_from_entropy
@@ -141,14 +154,14 @@ INLINE static float Til_soundspeed_from_entropy(float density, float entropy,
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// gas_entropy_from_internal_energy
INLINE static float Til_entropy_from_internal_energy(
float density, float u, const struct Til_params *mat) {
- return 0;
+ return 0.f;
}
// gas_pressure_from_internal_energy
@@ -156,35 +169,37 @@ 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 eta_sq = eta * eta;
const float mu = eta - 1.f;
const float nu = 1.f / eta - 1.f;
+ const float w = u / (mat->u_0 * eta_sq) + 1.f;
+ const float w_inv = 1.f / w;
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;
+ P_c = (mat->a + mat->b * w_inv) * 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)) *
+ (mat->b * density * u * w_inv + mat->A * mu * expf(-mat->beta * nu)) *
expf(-mat->alpha * nu * nu);
// Condensed or cold state
- if ((1.f < eta) || (u < mat->E_iv)) {
+ if ((1.f < eta) || (u < mat->u_iv)) {
P = P_c;
}
// Expanded and hot state
- else if ((eta < 1.f) && (mat->E_cv < u)) {
+ else if ((eta < 1.f) && (mat->u_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);
+ P = ((u - mat->u_iv) * P_e + (mat->u_cv - u) * P_c) /
+ (mat->u_cv - mat->u_iv);
}
// Minimum pressure
@@ -201,81 +216,78 @@ INLINE static float Til_internal_energy_from_pressure(
error("This EOS function is not yet implemented!");
- return 0;
+ return 0.f;
}
// 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;
+ const float rho_0_inv = 1.f / mat->rho_0;
+ const float eta = density * rho_0_inv;
+ const float rho_inv = 1.f / density;
+ const float eta_sq = eta * eta;
+ const float mu = eta - 1.f;
+ const float nu = 1.f / eta - 1.f;
+ const float w = u / (mat->u_0 * eta_sq) + 1.f;
+ const float w_inv = 1.f / w;
+ const float w_inv_sq = w_inv * w_inv;
+ const float exp_beta = expf(-mat->beta * nu);
+ const float exp_alpha = expf(-mat->alpha * nu * nu);
+ float P_c, P_e, c_sq_c, c_sq_e, c_sq;
+
+ // Condensed or cold
+ if (eta < mat->eta_min) {
+ P_c = 0.f;
+ } else {
+ P_c = (mat->a + mat->b * w_inv) * density * u + mat->A * mu +
+ mat->B * mu * mu;
+ }
+ c_sq_c = P_c * rho_inv * (1.f - mat->a - mat->b * w_inv) +
+ mat->b * (w - 1.f) * w_inv_sq * (2 * u + P_c * rho_inv) +
+ rho_inv * (mat->A + mat->B * (eta_sq - 1.f));
+
+ c_sq_c = fmax(c_sq_c, mat->A * rho_0_inv);
+
+ // Expanded and hot
+ P_e = mat->a * density * u +
+ (mat->b * density * u * w_inv + mat->A * mu * exp_beta) * exp_alpha;
+
+ c_sq_e = P_e * rho_inv * (1.f - mat->a) +
+ (mat->b * density * u / (w * w * eta_sq) *
+ (rho_inv / mat->u_0 * (2 * u - P_e * rho_inv * eta_sq) +
+ 2.f * mat->alpha * nu * rho_0_inv) +
+ mat->A * rho_0_inv *
+ (1 + mu / eta_sq * (mat->beta + 2.f * mat->alpha * nu - eta)) *
+ exp_beta) *
+ exp_alpha;
+
+ // Condensed or cold state
+ if ((1.f < eta) || (u < mat->u_iv)) {
+ c_sq = c_sq_c;
+ }
+ // Expanded and hot state
+ else if ((eta < 1.f) && (mat->u_cv < u)) {
+ c_sq = c_sq_e;
+ }
+ // Hybrid state
+ else {
+ c_sq = ((u - mat->u_iv) * c_sq_e + (mat->u_cv - u) * c_sq_c) /
+ (mat->u_cv - mat->u_iv);
+
+ c_sq = fmax(c_sq_c, mat->A * rho_0_inv);
+ }
+
+ return sqrtf(c_sq);
}
// 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);
+ error("This EOS function is not yet implemented!");
- return c;
+ return 0.f;
}
#endif /* SWIFT_TILLOTSON_EQUATION_OF_STATE_H */
diff --git a/src/hydro.h b/src/hydro.h
index 950f63526a1590fa0fdcf2bfb5e650a2dfe14431..b3716996cc4da68f9445adccd12315b32d81a34c 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -62,9 +62,9 @@
#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"
+#elif defined(PLANETARY_SPH)
+#include "./hydro/Planetary/hydro.h"
+#include "./hydro/Planetary/hydro_iact.h"
#define SPH_IMPLEMENTATION "Minimal version of SPH with multiple materials"
#else
#error "Invalid choice of SPH variant"
diff --git a/src/hydro/MinimalMultiMat/hydro.h b/src/hydro/Planetary/hydro.h
similarity index 92%
rename from src/hydro/MinimalMultiMat/hydro.h
rename to src/hydro/Planetary/hydro.h
index cfad6b2b2b389da9f423540cb30f1df4cebc5416..134893d9bb01e86de0a17e196d307e6972372c03 100644
--- a/src/hydro/MinimalMultiMat/hydro.h
+++ b/src/hydro/Planetary/hydro.h
@@ -17,11 +17,11 @@
* 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
+#ifndef SWIFT_PLANETARY_HYDRO_H
+#define SWIFT_PLANETARY_HYDRO_H
/**
- * @file MinimalMultiMat/hydro.h
+ * @file Planetary/hydro.h
* @brief Minimal conservative implementation of SPH (Non-neighbour loop
* equations) with multiple materials.
*
@@ -44,6 +44,12 @@
#include "kernel_hydro.h"
#include "minmax.h"
+/*
+ * Note: Define PLANETARY_SPH_BALSARA to use the Balsara (1995) switch for
+ * the artificial viscosity, instead of the default Monaghan (1992).
+ * i.e. compile with: make CFLAGS=-DPLANETARY_SPH_BALSARA to use.
+ */
+
/**
* @brief Returns the comoving internal energy of a particle
*
@@ -257,6 +263,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
p->u_dt = du_dt;
}
+
/**
* @brief Computes the hydro time-step of a given particle
*
@@ -391,23 +398,52 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {
+#ifdef PLANETARY_SPH_BALSARA
+ const float fac_mu = cosmo->a_factor_mu;
+
+ /* Compute the norm of the curl */
+ const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
+ p->density.rot_v[1] * p->density.rot_v[1] +
+ p->density.rot_v[2] * p->density.rot_v[2]);
+
+ /* Compute the norm of div v */
+ const float abs_div_v = fabsf(p->density.div_v);
+#endif // PLANETARY_SPH_BALSARA
+
/* 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);
+ gas_soundspeed_from_internal_energy(p->rho, p->u, 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);
+ float grad_h_term;
+ const float grad_h_term_inv =
+ 1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv;
+ /* Avoid 1/0 from only having one neighbour right at the edge of the kernel */
+ if (grad_h_term_inv != 0.f) {
+ grad_h_term = 1.f / grad_h_term_inv;
+ } else {
+ grad_h_term = 0.f;
+ }
+
+#ifdef PLANETARY_SPH_BALSARA
+ /* Compute the Balsara switch */
+ const float balsara =
+ abs_div_v / (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
+#endif // PLANETARY_SPH_BALSARA
/* Update variables. */
p->force.f = grad_h_term;
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+#ifdef PLANETARY_SPH_BALSARA
+ p->force.balsara = balsara;
+#endif // PLANETARY_SPH_BALSARA
}
/**
@@ -494,7 +530,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Compute the new sound speed */
const float soundspeed =
- gas_soundspeed_from_pressure(p->rho, pressure, p->mat_id);
+ gas_soundspeed_from_internal_energy(p->rho, p->u, p->mat_id);
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
@@ -631,4 +667,4 @@ hydro_set_init_internal_energy(struct part *p, float u_init) {
p->u = u_init;
}
-#endif /* SWIFT_MINIMAL_MULTI_MAT_HYDRO_H */
+#endif /* SWIFT_PLANETARY_HYDRO_H */
diff --git a/src/hydro/MinimalMultiMat/hydro_debug.h b/src/hydro/Planetary/hydro_debug.h
similarity index 89%
rename from src/hydro/MinimalMultiMat/hydro_debug.h
rename to src/hydro/Planetary/hydro_debug.h
index 17b624ad0f660152be4ba685905a3c855e1761f8..74261f3b49e2881af1c403013005560efa53a7f1 100644
--- a/src/hydro/MinimalMultiMat/hydro_debug.h
+++ b/src/hydro/Planetary/hydro_debug.h
@@ -17,13 +17,12 @@
* 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
+#ifndef SWIFT_PLANETARY_HYDRO_DEBUG_H
+#define SWIFT_PLANETARY_HYDRO_DEBUG_H
/**
- * @file MinimalMultiMat/hydro_debug.h
- * @brief MinimalMultiMat conservative implementation of SPH (Debugging
- * routines)
+ * @file Planetary/hydro_debug.h
+ * @brief Minimal 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
@@ -51,4 +50,4 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
p->density.wcount, p->rho, p->density.rho_dh, p->time_bin, p->mat_id);
}
-#endif /* SWIFT_MINIMAL_MULTI_MAT_HYDRO_DEBUG_H */
+#endif /* SWIFT_PLANETARY_HYDRO_DEBUG_H */
diff --git a/src/hydro/MinimalMultiMat/hydro_iact.h b/src/hydro/Planetary/hydro_iact.h
similarity index 89%
rename from src/hydro/MinimalMultiMat/hydro_iact.h
rename to src/hydro/Planetary/hydro_iact.h
index 5984c1c483546d87800792ced0ffcc41e0aaa408..bf96034696806e3adff1d8ba7f385af65461b9ea 100644
--- a/src/hydro/MinimalMultiMat/hydro_iact.h
+++ b/src/hydro/Planetary/hydro_iact.h
@@ -17,13 +17,12 @@
* 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
+#ifndef SWIFT_PLANETARY_HYDRO_IACT_H
+#define SWIFT_PLANETARY_HYDRO_IACT_H
/**
- * @file MinimalMultiMat/hydro_iact.h
- * @brief MinimalMultiMat conservative implementation of SPH (Neighbour loop
- * equations)
+ * @file Planetary/hydro_iact.h
+ * @brief Minimal 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
@@ -177,7 +176,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
(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 ? */
+#ifdef PLANETARY_SPH_BALSARA
+ /* Balsara term */
+ const float balsara_i = pi->force.balsara;
+ const float balsara_j = pj->force.balsara;
+#endif // PLANETARY_SPH_BALSARA
+
+ /* Are the particles moving towards each other? */
const float omega_ij = min(dvdr, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
@@ -186,9 +191,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float cj = pj->force.soundspeed;
const float v_sig = ci + cj - 3.f * mu_ij;
- /* Construct the full viscosity term */
+ /* Now construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
+#ifdef PLANETARY_SPH_BALSARA
+ const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
+ (balsara_i + balsara_j) / rho_ij;
+#else
const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+#endif // PLANETARY_SPH_BALSARA
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
@@ -220,7 +230,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
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 */
+ /* Internal energy time derivative */
pi->u_dt += du_dt_i * mj;
pj->u_dt += du_dt_j * mi;
@@ -290,18 +300,31 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
(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 ? */
+#ifdef PLANETARY_SPH_BALSARA
+ /* Balsara term */
+ const float balsara_i = pi->force.balsara;
+ const float balsara_j = pj->force.balsara;
+#endif // PLANETARY_SPH_BALSARA
+
+ /* Are the particles moving towards each other? */
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 */
+ /* Compute sound speeds */
const float ci = pi->force.soundspeed;
const float cj = pj->force.soundspeed;
+
+ /* Signal velocity */
const float v_sig = ci + cj - 3.f * mu_ij;
/* Construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
+#ifdef PLANETARY_SPH_BALSARA
+ const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
+ (balsara_i + balsara_j) / rho_ij;
+#else
const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+#endif // PLANETARY_SPH_BALSARA
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
@@ -327,7 +350,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term;
- /* Internal energy time derivatibe */
+ /* Internal energy time derivative */
pi->u_dt += du_dt_i * mj;
/* Get the time derivative for h. */
@@ -337,4 +360,4 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.v_sig = max(pi->force.v_sig, v_sig);
}
-#endif /* SWIFT_MINIMAL_MULTI_MAT_HYDRO_IACT_H */
+#endif /* SWIFT_PLANETARY_HYDRO_IACT_H */
diff --git a/src/hydro/MinimalMultiMat/hydro_io.h b/src/hydro/Planetary/hydro_io.h
similarity index 95%
rename from src/hydro/MinimalMultiMat/hydro_io.h
rename to src/hydro/Planetary/hydro_io.h
index 7f41f5e227b6c8a8904b5546a2568b4700109abd..afb37d884494fd02e30c143194804a2b49a77be0 100644
--- a/src/hydro/MinimalMultiMat/hydro_io.h
+++ b/src/hydro/Planetary/hydro_io.h
@@ -17,12 +17,12 @@
* 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
+#ifndef SWIFT_PLANETARY_HYDRO_IO_H
+#define SWIFT_PLANETARY_HYDRO_IO_H
/**
- * @file MinimalMultiMat/hydro_io.h
- * @brief MinimalMultiMat conservative implementation of SPH (i/o routines)
+ * @file Planetary/hydro_io.h
+ * @brief Minimal 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
@@ -197,8 +197,14 @@ INLINE static 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");
+#ifdef PLANETARY_SPH_BALSARA
+ io_write_attribute_s(
+ h_grpsph, "Viscosity Model",
+ "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
+#else
io_write_attribute_s(h_grpsph, "Viscosity Model",
"Minimal treatment as in Monaghan (1992)");
+#endif // PLANETARY_SPH_BALSARA
/* Time integration properties */
io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
@@ -212,4 +218,4 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) {
*/
INLINE static int writeEntropyFlag(void) { return 0; }
-#endif /* SWIFT_MINIMAL_MULTI_MAT_HYDRO_IO_H */
+#endif /* SWIFT_PLANETARY_HYDRO_IO_H */
diff --git a/src/hydro/MinimalMultiMat/hydro_part.h b/src/hydro/Planetary/hydro_part.h
similarity index 90%
rename from src/hydro/MinimalMultiMat/hydro_part.h
rename to src/hydro/Planetary/hydro_part.h
index dad13e889aa531636e34846825109086177b3dae..7d1fc8f6729992bfdf2eeaba6e33cc9a7b071655 100644
--- a/src/hydro/MinimalMultiMat/hydro_part.h
+++ b/src/hydro/Planetary/hydro_part.h
@@ -17,13 +17,12 @@
* 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
+#ifndef SWIFT_PLANETARY_HYDRO_PART_H
+#define SWIFT_PLANETARY_HYDRO_PART_H
/**
- * @file MinimalMultiMat/hydro_part.h
- * @brief MinimalMultiMat conservative implementation of SPH (Particle
- * definition)
+ * @file Planetary/hydro_part.h
+ * @brief Minimal 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
@@ -127,6 +126,14 @@ struct part {
/*! Derivative of density with respect to h */
float rho_dh;
+#ifdef PLANETARY_SPH_BALSARA
+ /*! Velocity divergence. */
+ float div_v;
+
+ /*! Velocity curl. */
+ float rot_v[3];
+#endif // PLANETARY_SPH_BALSARA
+
} density;
/**
@@ -153,6 +160,11 @@ struct part {
/*! Time derivative of smoothing length */
float h_dt;
+#ifdef PLANETARY_SPH_BALSARA
+ /*! Balsara switch */
+ float balsara;
+#endif // PLANETARY_SPH_BALSARA
+
} force;
};
@@ -177,4 +189,4 @@ struct part {
} SWIFT_STRUCT_ALIGN;
-#endif /* SWIFT_MINIMAL_MULTI_MAT_HYDRO_PART_H */
+#endif /* SWIFT_PLANETARY_HYDRO_PART_H */
diff --git a/src/hydro_io.h b/src/hydro_io.h
index d752bb8bc03f619fe759fc8f5de32a01b3a61abe..b6e0c36cc7415a1f628a109795aa98b4f583036c 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -39,8 +39,8 @@
#include "./hydro/GizmoMFM/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"
+#elif defined(PLANETARY_SPH)
+#include "./hydro/Planetary/hydro_io.h"
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index f79fd832248fba8fbc55bd9fcec57e645be93159..905bf6973447b1ddd1c174b2e65d6841917ef736 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -135,6 +135,14 @@ void hydro_props_init(struct hydro_props *p,
mean_molecular_weight = 4. / (1. + 3. * p->hydrogen_mass_fraction);
p->minimal_internal_energy = u_min / mean_molecular_weight;
+
+#ifdef PLANETARY_SPH
+#ifdef PLANETARY_SPH_BALSARA
+ message("Planetary SPH: Balsara switch enabled");
+#else
+ message("Planetary SPH: Balsara switch disabled");
+#endif // PLANETARY_SPH_BALSARA
+#endif // PLANETARY_SPH
}
/**
diff --git a/src/multipole.h b/src/multipole.h
index e0e6da32a2950d7fce164b2abc422302b7c7de5e..b9407d7698146520fed8331697fd539df5a0780a 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -1037,6 +1037,13 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
vel[2] += gparts[k].v_full[2] * m;
}
+#ifdef PLANETARY_SPH
+ /* Prevent FPE from zero mass with the temporary outside-the-box particles */
+ if (mass == 0.f) {
+ mass = FLT_MIN;
+ }
+#endif // PLANETARY_SPH
+
/* Final operation on CoM */
const double imass = 1.0 / mass;
com[0] *= imass;
diff --git a/src/outputlist.c b/src/outputlist.c
index 782bdeb3eb53aeb1c259ca0283c8ccaa15d68949..fd33370ca45f25c17ecd2cc8df622138842507f3 100644
--- a/src/outputlist.c
+++ b/src/outputlist.c
@@ -208,7 +208,8 @@ void output_list_read_next_time(struct output_list *t, const struct engine *e,
* time)
*/
void output_list_init(struct output_list **list, const struct engine *e,
- const char *name, double *delta_time, double *time_first) {
+ const char *name, double *delta_time,
+ double *time_first) {
struct swift_params *params = e->parameter_file;
/* get cosmo */
diff --git a/src/parallel_io.c b/src/parallel_io.c
index b82443b33ba767cfe2050cf300535db924bfb537..a3a71d0ac5703f194bfa8e53c9b0b62b25e6a602 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -928,10 +928,11 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
char fileName[FILENAME_BUFFER_SIZE];
if (e->snapshot_label_delta == 1)
snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
- e->snapshot_output_count);
+ e->snapshot_output_count + e->snapshot_label_first);
else
snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%06i.hdf5", baseName,
- e->snapshot_output_count * e->snapshot_label_delta);
+ e->snapshot_output_count * e->snapshot_label_delta +
+ e->snapshot_label_first);
/* Open HDF5 file with the chosen parameters */
hid_t h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
@@ -1204,8 +1205,13 @@ void write_output_parallel(struct engine* e, const char* baseName,
/* HDF5 File name */
char fileName[FILENAME_BUFFER_SIZE];
- snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
- e->snapshot_output_count);
+ if (e->snapshot_label_delta == 1)
+ snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
+ e->snapshot_output_count + e->snapshot_label_first);
+ else
+ snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%06i.hdf5", baseName,
+ e->snapshot_output_count * e->snapshot_label_delta +
+ e->snapshot_label_first);
/* Prepare some file-access properties */
hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
diff --git a/src/part.h b/src/part.h
index bca84cc0212e79e15ffbeeeb0bbcfc714d5481be..145bf2111771d8ad254affb213b93b7ac829f1a6 100644
--- a/src/part.h
+++ b/src/part.h
@@ -69,8 +69,8 @@
#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"
+#elif defined(PLANETARY_SPH)
+#include "./hydro/Planetary/hydro_part.h"
#define hydro_need_extra_init_loop 0
#else
#error "Invalid choice of SPH variant"
diff --git a/src/runner.c b/src/runner.c
index 7771e247a07866297576856e4f5099fa8a1e55da..57cfe863d2a1efbe8c62c91d2d4b02d495630c22 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -766,6 +766,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Double h and try again */
h_new = 2.f * h_old;
+
} else {
/* Finish the density calculation */
@@ -780,6 +781,46 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
p->density.wcount_dh * h_old_dim +
hydro_dimension * p->density.wcount * h_old_dim_minus_one;
+ /* Skip if h is already h_max and we don't have enough neighbours */
+ if ((p->h >= hydro_h_max) && (f < 0.f)) {
+
+ /* We have a particle whose smoothing length is already set (wants to
+ * be larger but has already hit the maximum). So, just tidy up as if
+ * the smoothing length had converged correctly */
+
+#ifdef EXTRA_HYDRO_LOOP
+
+ /* As of here, particle gradient variables will be set. */
+ /* The force variables are set in the extra ghost. */
+
+ /* Compute variables required for the gradient loop */
+ hydro_prepare_gradient(p, xp, cosmo);
+
+ /* The particle gradient values are now set. Do _NOT_
+ try to read any particle density variables! */
+
+ /* Prepare the particle for the gradient loop over neighbours */
+ hydro_reset_gradient(p);
+
+#else
+ /* As of here, particle force variables will be set. */
+
+ /* Compute variables required for the force loop */
+ hydro_prepare_force(p, xp, cosmo);
+
+ /* The particle force values are now set. Do _NOT_
+ try to read any particle density variables! */
+
+ /* Prepare the particle for the force loop over neighbours */
+ hydro_reset_acceleration(p);
+
+#endif /* EXTRA_HYDRO_LOOP */
+
+ continue;
+ }
+
+ /* Normal case: Use Newton-Raphson to get a better value of h */
+
/* Avoid floating point exception from f_prime = 0 */
h_new = h_old - f / (f_prime + FLT_MIN);
diff --git a/src/serial_io.c b/src/serial_io.c
index dafa75ab0baacb1b5ddeee34020c9773893bced7..3a7d2e5a68873ca9523fe09bbf19fb2e185482dd 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -759,10 +759,11 @@ void write_output_serial(struct engine* e, const char* baseName,
char fileName[FILENAME_BUFFER_SIZE];
if (e->snapshot_label_delta == 1)
snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
- e->snapshot_output_count);
+ e->snapshot_output_count + e->snapshot_label_first);
else
snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%06i.hdf5", baseName,
- e->snapshot_output_count * e->snapshot_label_delta);
+ e->snapshot_output_count * e->snapshot_label_delta +
+ e->snapshot_label_first);
/* Compute offset in the file and total number of particles */
size_t N[swift_type_count] = {Ngas, Ndm, 0, 0, Nstars, 0};
diff --git a/src/single_io.c b/src/single_io.c
index a0f02878b52c89beca94d15c09ef7d456ce0a4eb..8cbb6743d38a022581273a0a0b03c9b3b6fda32e 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -625,10 +625,11 @@ void write_output_single(struct engine* e, const char* baseName,
char fileName[FILENAME_BUFFER_SIZE];
if (e->snapshot_label_delta == 1)
snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
- e->snapshot_output_count);
+ e->snapshot_output_count + e->snapshot_label_first);
else
snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%06i.hdf5", baseName,
- e->snapshot_output_count * e->snapshot_label_delta);
+ e->snapshot_output_count * e->snapshot_label_delta +
+ e->snapshot_label_first);
/* First time, we need to create the XMF file */
if (e->snapshot_output_count == 0) xmf_create_file(baseName);
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 70af9dfb19aeee812fdc90733a604f795b7f478f..2a2c20dbb064539b481e169b49b74389e79a8174 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -120,7 +120,7 @@ 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) || defined(HOPKINS_PU_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
-#elif defined(MINIMAL_MULTI_MAT_SPH)
+#elif defined(PLANETARY_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
part->primitives.P = pressure;
@@ -407,8 +407,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(MINIMAL_MULTI_MAT_SPH) || \
- defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH) || \
+#if defined(MINIMAL_SPH) || defined(PLANETARY_SPH) || \
+ defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH) || \
defined(HOPKINS_PU_SPH)
0.f,
#else
diff --git a/tests/testEOS.c b/tests/testEOS.c
index d090d83d77a16b2c5b4506c5f9224b3e4434d1be..4a1e666b47acc55a5ed7f1800e7199a1abb5e821 100644
--- a/tests/testEOS.c
+++ b/tests/testEOS.c
@@ -74,6 +74,11 @@
* P_1_0 ... ... P_1_num_u
* ... ... ... ...
* P_num_rho_0 ... P_num_rho_num_u
+ * c_0_0 c_0_1 ... c_0_num_u # Array of sound speeds, c(rho,
+ * u)
+ * c_1_0 ... ... c_1_num_u
+ * ... ... ... ...
+ * c_num_rho_0 ... c_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
@@ -83,21 +88,24 @@
#ifdef EOS_PLANETARY
int main(int argc, char *argv[]) {
- float rho, log_rho, log_u, P;
+ float rho, u, log_rho, log_u, P, c;
struct unit_system us;
+ struct swift_params *params =
+ (struct swift_params *)malloc(sizeof(struct swift_params));
+ if (params == NULL) error("Error allocating memory for the parameter file.");
const struct phys_const *phys_const = 0; // Unused placeholder
- 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),
+ float log_rho_min = logf(1e-4f), log_rho_max = logf(1e3f), // Densities (cgs)
+ log_u_min = logf(1e4f),
+ log_u_max = logf(1e13f), // Sp. int. energies (SI)
+ 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;
+ int mat_id_in, do_output;
// Default sys args
const int mat_id_def = eos_planetary_id_HM80_ice;
const int do_output_def = 0;
@@ -106,34 +114,40 @@ int main(int argc, char *argv[]) {
switch (argc) {
case 1:
// Default both
- mat_id = mat_id_def;
+ mat_id_in = mat_id_def;
do_output = do_output_def;
break;
case 2:
// Read mat_id, default do_output
- mat_id = atoi(argv[1]);
+ mat_id_in = atoi(argv[1]);
do_output = do_output_def;
break;
case 3:
// Read both
- mat_id = atoi(argv[1]);
+ mat_id_in = 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
+ mat_id_in = mat_id_def; // Ignored, just here to keep the compiler happy
do_output = do_output_def;
};
+ enum eos_planetary_material_id mat_id =
+ (enum eos_planetary_material_id)mat_id_in;
+
/* Greeting message */
printf("This is %s\n", package_description());
// Check material ID
- // Material base type
- switch ((int)(mat_id / eos_planetary_type_factor)) {
+ const enum eos_planetary_type_id type =
+ (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
+
+ // Select the material base type
+ switch (type) {
// Tillotson
case eos_planetary_type_Til:
switch (mat_id) {
@@ -174,27 +188,23 @@ int main(int argc, char *argv[]) {
};
break;
- // ANEOS
- case eos_planetary_type_ANEOS:
+ // SESAME
+ case eos_planetary_type_SESAME:
switch (mat_id) {
- case eos_planetary_id_ANEOS_iron:
- printf(" ANEOS iron \n");
+ case eos_planetary_id_SESAME_iron:
+ printf(" SESAME basalt 7530 \n");
break;
- case eos_planetary_id_MANEOS_forsterite:
- printf(" MANEOS forsterite \n");
+ case eos_planetary_id_SESAME_basalt:
+ printf(" SESAME basalt 7530 \n");
break;
- default:
- error("Unknown material ID! mat_id = %d \n", mat_id);
- };
- break;
+ case eos_planetary_id_SESAME_water:
+ printf(" SESAME water 7154 \n");
+ break;
- // SESAME
- case eos_planetary_type_SESAME:
- switch (mat_id) {
- case eos_planetary_id_SESAME_iron:
- printf(" SESAME iron \n");
+ case eos_planetary_id_SS08_water:
+ printf(" Senft & Stewart (2008) SESAME-like water \n");
break;
default:
@@ -206,8 +216,11 @@ int main(int argc, char *argv[]) {
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);
+ // Convert to internal units
+ // Earth masses and radii
+ // units_init(&us, 5.9724e27, 6.3710e8, 1.f, 1.f, 1.f);
+ // SI
+ units_init(&us, 1000.f, 100.f, 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(
@@ -215,11 +228,51 @@ int main(int argc, char *argv[]) {
log_u_max += logf(J_kg_to_erg_g / units_cgs_conversion_factor(
&us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
+ // Set the input parameters
+ // Which EOS to initialise
+ parser_set_param(params, "EoS:planetary_use_Til:1");
+ parser_set_param(params, "EoS:planetary_use_HM80:1");
+ parser_set_param(params, "EoS:planetary_use_SESAME:1");
+ // Table file names
+ parser_set_param(params,
+ "EoS:planetary_HM80_HHe_table_file:"
+ "../examples/planetary_HM80_HHe.txt");
+ parser_set_param(params,
+ "EoS:planetary_HM80_ice_table_file:"
+ "../examples/planetary_HM80_ice.txt");
+ parser_set_param(params,
+ "EoS:planetary_HM80_rock_table_file:"
+ "../examples/planetary_HM80_rock.txt");
+ parser_set_param(params,
+ "EoS:planetary_SESAME_iron_table_file:"
+ "../examples/planetary_SESAME_iron_2140.txt");
+ parser_set_param(params,
+ "EoS:planetary_SESAME_basalt_table_file:"
+ "../examples/planetary_SESAME_basalt_7530.txt");
+ parser_set_param(params,
+ "EoS:planetary_SESAME_water_table_file:"
+ "../examples/planetary_SESAME_water_7154.txt");
+ parser_set_param(params,
+ "EoS:planetary_SS08_water_table_file:"
+ "../examples/planetary_SS08_water.txt");
+
// Initialise the EOS materials
eos_init(&eos, phys_const, &us, params);
+ // Manual debug testing
+ if (1) {
+ printf("\n ### MANUAL DEBUG TESTING ### \n");
+
+ rho = 5960;
+ u = 1.7e8;
+ P = gas_pressure_from_internal_energy(rho, u, eos_planetary_id_HM80_ice);
+ printf("u = %.2e, rho = %.2e, P = %.2e \n", u, rho, P);
+
+ return 0;
+ }
+
// Output file
- sprintf(filename, "testEOS_rho_u_P_%d.txt", mat_id);
+ sprintf(filename, "testEOS_rho_u_P_c_%d.txt", mat_id);
FILE *f = fopen(filename, "w");
if (f == NULL) {
printf("Could not open output file!\n");
@@ -270,6 +323,21 @@ int main(int argc, char *argv[]) {
if (do_output == 1) fprintf(f, "\n");
}
+
+ // Sound speeds
+ for (int i = 0; i < num_rho; i++) {
+ rho = A1_rho[i];
+
+ for (int j = 0; j < num_u; j++) {
+ c = gas_soundspeed_from_internal_energy(rho, A1_u[j], mat_id);
+
+ if (do_output == 1)
+ fprintf(f, "%.6e ",
+ c * units_cgs_conversion_factor(&us, UNIT_CONV_SPEED));
+ }
+
+ if (do_output == 1) fprintf(f, "\n");
+ }
fclose(f);
return 0;
diff --git a/tests/testEOS.py b/tests/testEOS.py
index 363bab200b58c65fa24cc033af4b8d3c04b7b503..a2a31a248a2073a834d9543b706a6a12ba12796c 100644
--- a/tests/testEOS.py
+++ b/tests/testEOS.py
@@ -18,8 +18,8 @@
#
##############################################################################
"""
-Plot the output of testEOS to show how the equation of state pressure varies
-with density and specific internal energy.
+Plot the output of testEOS to show how the equation of state pressure and sound
+speed varies with density and specific internal energy.
Usage:
python testEOS.py (mat_id)
@@ -37,9 +37,13 @@ Text file contains:
P_1_0 ... ... P_1_num_u
... ... ... ...
P_num_rho_0 ... P_num_rho_num_u
+ c_0_0 c_0_1 ... c_0_num_u # Array of sound speeds, c(rho, u)
+ c_1_0 ... ... c_1_num_u
+ ... ... ... ...
+ c_num_rho_0 ... c_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
+designed for (e.g. outside table limits), to help test the EOS in cases of
unexpected particle behaviour.
"""
@@ -58,8 +62,7 @@ type_factor = 100
Di_type = {
'Til' : 1,
'HM80' : 2,
- 'ANEOS' : 3,
- 'SESAME' : 4,
+ 'SESAME' : 3,
}
Di_material = {
# Tillotson
@@ -70,11 +73,11 @@ Di_material = {
'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,
+ 'SESAME_iron' : Di_type['SESAME']*type_factor, # 2140
+ 'SESAME_basalt' : Di_type['SESAME']*type_factor + 1, # 7530
+ 'SESAME_water' : Di_type['SESAME']*type_factor + 2, # 7154
+ 'SS08_water' : Di_type['SESAME']*type_factor + 3, # Senft & Stewart (2008)
}
# Invert so the mat_id are the keys
Di_mat_id = {mat_id : mat for mat, mat_id in Di_material.iteritems()}
@@ -82,6 +85,7 @@ 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
+cm_s_to_m_s = 1e-2
if __name__ == '__main__':
# Sys args
@@ -101,7 +105,7 @@ if __name__ == '__main__':
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
+ filename = "testEOS_rho_u_P_c_%d.txt" % mat_id
# Load the header info and density and energy arrays
with open(filename) as f:
@@ -110,31 +114,37 @@ if __name__ == '__main__':
A1_rho = np.array(f.readline().split(), dtype=float)
A1_u = np.array(f.readline().split(), dtype=float)
- # Load pressure array
+ # Load pressure and soundspeed arrays
A2_P = np.loadtxt(filename, skiprows=4)
+ A2_c = A2_P[num_rho:]
+ A2_P = A2_P[:num_rho]
- # 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
+ # Convert pressures from cgs (Barye) to Mbar
+ A2_P *= Ba_to_Mbar
+ # Convert sound speeds from cgs to SI
+ A1_u *= cm_s_to_m_s
# 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)
+ assert A2_c.shape == (num_rho, num_u)
# Plot
+ # Pressure: P(rho) at fixed u
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))
+ A1_colour = matplotlib.cm.rainbow(np.linspace(0, 1, num_u))
- for i, idx in enumerate(A1_idx):
- plt.plot(A1_rho, A2_P[:, idx], c=A1_colour[i],
- label=r"%.2e" % A1_u[idx])
+ for i_u, u in enumerate(A1_u):
+ if i_u%10 == 0:
+ plt.plot(A1_rho, A2_P[:, i_u], c=A1_colour[i_u],
+ label=r"%.2e" % u)
+ else:
+ plt.plot(A1_rho, A2_P[:, i_u], c=A1_colour[i_u])
plt.legend(title="Sp. Int. Energy (J kg$^{-1}$)")
plt.xscale('log')
@@ -144,7 +154,31 @@ if __name__ == '__main__':
plt.title(mat)
plt.tight_layout()
- plt.savefig("testEOS_%d.png" % mat_id)
+ plt.savefig("testEOS_P_%d.png" % mat_id)
+ plt.close()
+
+ # Sound speed: c(rho) at fixed u
+ plt.figure(figsize=(7, 7))
+ ax = plt.gca()
+
+ A1_colour = matplotlib.cm.rainbow(np.linspace(0, 1, num_u))
+
+ for i_u, u in enumerate(A1_u):
+ if i_u%10 == 0:
+ plt.plot(A1_rho, A2_c[:, i_u], c=A1_colour[i_u],
+ label=r"%.2e" % u)
+ else:
+ plt.plot(A1_rho, A2_c[:, i_u], c=A1_colour[i_u])
+
+ 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"Sound Speed (m s^{-1})")
+ plt.title(mat)
+ plt.tight_layout()
+
+ plt.savefig("testEOS_c_%d.png" % mat_id)
plt.close()
diff --git a/tests/testEOS.sh b/tests/testEOS.sh
index 411ac746be186bfe5758e03c2a852e081daefd10..bcd87eabbf15a962808843dda76d1829f2917c97 100755
--- a/tests/testEOS.sh
+++ b/tests/testEOS.sh
@@ -13,6 +13,10 @@ A1_mat_id=(
200
201
202
+ 300
+ 301
+ 302
+ 303
)
for mat_id in "${A1_mat_id[@]}"
diff --git a/tests/testEOS_plot.sh b/tests/testEOS_plot.sh
index 39108c5e19d8f4474de508e205951a1fd0aebcc9..5fd7f4976496223e467aae65b2846a8c4e1b7485 100755
--- a/tests/testEOS_plot.sh
+++ b/tests/testEOS_plot.sh
@@ -2,6 +2,8 @@
echo ""
+rm -f testEOS*.png
+
echo "Plotting testEOS output for each planetary material"
A1_mat_id=(
@@ -11,6 +13,10 @@ A1_mat_id=(
200
201
202
+ 300
+ 301
+ 302
+ 303
)
for mat_id in "${A1_mat_id[@]}"
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index 0a7354f0d2a5e1853ba2c22d696dbb910de2b667..306f14a35ca047430f67e33e9fd63848e9207b68 100644
--- a/tests/testInteractions.c
+++ b/tests/testInteractions.c
@@ -107,8 +107,8 @@ struct part *make_particles(size_t count, double *offset, double spacing,
*/
void prepare_force(struct part *parts, size_t count) {
-#if !defined(GIZMO_MFV_SPH) && !defined(SHADOWFAX_SPH) && \
- !defined(MINIMAL_SPH) && !defined(MINIMAL_MULTI_MAT_SPH) && \
+#if !defined(GIZMO_MFV_SPH) && !defined(SHADOWFAX_SPH) && \
+ !defined(MINIMAL_SPH) && !defined(PLANETARY_SPH) && \
!defined(HOPKINS_PU_SPH)
struct part *p;
for (size_t i = 0; i < count; ++i) {
@@ -136,8 +136,7 @@ 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(MINIMAL_MULTI_MAT_SPH) || \
- defined(SHADOWFAX_SPH)
+#if defined(MINIMAL_SPH) || defined(PLANETARY_SPH) || defined(SHADOWFAX_SPH)
0.f,
#else
p->density.div_v,