diff --git a/.gitignore b/.gitignore
index 6b03972ab8a13017ae80e284771618d4c65f6c24..b7f8ac3a1466459891366cd774e922833b1cca74 100644
--- a/.gitignore
+++ b/.gitignore
@@ -114,6 +114,7 @@ tests/testPotentialSelf
tests/testPotentialPair
tests/testEOS
tests/testEOS*.txt
+tests/testEOS*.png
theory/latex/swift.pdf
theory/SPH/Kernels/kernels.pdf
diff --git a/tests/Makefile.am b/tests/Makefile.am
index ded168984750e1c9b1ecb0d5f94bf4094c283f2d..5789b1baf64554d5c53aa6696ad5727bcfda07cb 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -27,7 +27,7 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
testMatrixInversion testThreadpool testDump testLogger testInteractions.sh \
testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
- testPotentialPair
+ testPotentialPair testEOS
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
@@ -37,7 +37,7 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testAdiabaticIndex testRiemannExact testRiemannTRRS \
testRiemannHLLC testMatrixInversion testDump testLogger \
testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
- testGravityDerivatives testPotentialSelf testPotentialPair
+ testGravityDerivatives testPotentialSelf testPotentialPair testEOS
# Rebuild tests when SWIFT is updated.
$(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -105,6 +105,8 @@ testPotentialSelf_SOURCES = testPotentialSelf.c
testPotentialPair_SOURCES = testPotentialPair.c
+testEOS_SOURCES = testEOS.c
+
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testActivePair.sh \
test27cells.sh test27cellsPerturbed.sh testParser.sh testPeriodicBC.sh \
@@ -112,4 +114,5 @@ EXTRA_DIST = testReading.sh makeInput.py testActivePair.sh \
difffloat.py tolerance_125_normal.dat tolerance_125_perturbed.dat \
tolerance_27_normal.dat tolerance_27_perturbed.dat tolerance_27_perturbed_h.dat tolerance_27_perturbed_h2.dat \
tolerance_testInteractions.dat tolerance_pair_active.dat tolerance_pair_force_active.dat \
- fft_params.yml tolerance_periodic_BC_normal.dat tolerance_periodic_BC_perturbed.dat
+ fft_params.yml tolerance_periodic_BC_normal.dat tolerance_periodic_BC_perturbed.dat \
+ testEOS.sh testEOS_plot.sh
diff --git a/tests/testEOS.c b/tests/testEOS.c
new file mode 100644
index 0000000000000000000000000000000000000000..3f8ca496507268f2abf6875435b0042e7337adc9
--- /dev/null
+++ b/tests/testEOS.c
@@ -0,0 +1,282 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2018 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <fenv.h>
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+#include <pthread.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+
+/* Conditional headers. */
+#ifdef HAVE_LIBZ
+#include <zlib.h>
+#endif
+
+/* Local headers. */
+#include "equation_of_state.h"
+#include "swift.h"
+
+/* Engine policy flags. */
+#ifndef ENGINE_POLICY
+#define ENGINE_POLICY engine_policy_none
+#endif
+
+/**
+ * @brief Write a list of densities, energies, and resulting pressures to file
+ * from an equation of state.
+ *
+ * WORK IN PROGRESS
+ *
+ * So far only does the Hubbard & MacFarlane (1980) equations of state.
+ *
+ * Usage:
+ * $ ./testEOS (mat_id) (do_output)
+ *
+ * Sys args (optional):
+ * mat_id | int | Material ID, see equation_of_state.h for the options.
+ * Default: 201 (= id_HM80_ice).
+ *
+ * do_output | int | Set 1 to write the output file of rho, u, P values,
+ * set 0 for no output. Default: 0.
+ *
+ * Output text file contains:
+ * header
+ * num_rho num_u mat_id # Header values
+ * rho_0 rho_1 rho_2 ... rho_num_rho # Array of densities, rho
+ * u_0 u_1 u_2 ... u_num_u # Array of energies, u
+ * P_0_0 P_0_1 ... P_0_num_u # Array of pressures, P(rho, u)
+ * P_1_0 ... ... P_1_num_u
+ * ... ... ... ...
+ * P_num_rho_0 ... P_num_rho_num_u
+ *
+ * Note that the values tested extend beyond the range that most EOS are
+ * designed for (e.g. outside table limits), to help test the EOS in case of
+ * unexpected particle behaviour.
+ *
+ */
+
+#ifdef EOS_PLANETARY
+int main(int argc, char *argv[]) {
+ float rho, log_rho, log_u, P;
+ struct unit_system us;
+ const struct phys_const *phys_const = 0; // Unused placeholder
+ const struct swift_params *params = 0; // Unused placeholder
+ const float J_kg_to_erg_g = 1e4; // Convert J/kg to erg/g
+ char filename[64];
+ // Output table params
+ const int num_rho = 100, num_u = 100;
+ float log_rho_min = logf(1e-4), log_rho_max = logf(30.f),
+ log_u_min = logf(1e4), log_u_max = logf(1e10),
+ log_rho_step = (log_rho_max - log_rho_min) / (num_rho - 1.f),
+ log_u_step = (log_u_max - log_u_min) / (num_u - 1.f);
+ float A1_rho[num_rho], A1_u[num_u];
+ // Sys args
+ int mat_id, do_output;
+ // Default sys args
+ const int mat_id_def = id_HM80_ice;
+ const int do_output_def = 0;
+
+ // Check the number of system arguments and use defaults if not provided
+ switch (argc) {
+ case 1:
+ // Default both
+ mat_id = mat_id_def;
+ do_output = do_output_def;
+ break;
+
+ case 2:
+ // Read mat_id, default do_output
+ mat_id = atoi(argv[1]);
+ do_output = do_output_def;
+ break;
+
+ case 3:
+ // Read both
+ mat_id = atoi(argv[1]);
+ do_output = atoi(argv[2]);
+ break;
+
+ default:
+ error("Invalid number of system arguments!\n");
+ mat_id = mat_id_def; // Ignored, just here to keep the compiler happy
+ do_output = do_output_def;
+ };
+
+ /* Greeting message */
+ printf("This is %s\n", package_description());
+
+ // Check material ID
+ // Material base type
+ switch ((int)(mat_id / type_factor)) {
+ // Tillotson
+ case type_Til:
+ switch (mat_id) {
+ case id_Til_iron:
+ printf(" Tillotson iron \n");
+ break;
+
+ case id_Til_granite:
+ printf(" Tillotson granite \n");
+ break;
+
+ case id_Til_water:
+ printf(" Tillotson water \n");
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d \n", mat_id);
+ };
+ break;
+
+ // Hubbard & MacFarlane (1980)
+ case type_HM80:
+ switch (mat_id) {
+ case id_HM80_HHe:
+ printf(" Hubbard & MacFarlane (1980) hydrogen-helium atmosphere \n");
+ break;
+
+ case id_HM80_ice:
+ printf(" Hubbard & MacFarlane (1980) ice mix \n");
+ break;
+
+ case id_HM80_rock:
+ printf(" Hubbard & MacFarlane (1980) rock mix \n");
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d \n", mat_id);
+ };
+ break;
+
+ // ANEOS
+ case type_ANEOS:
+ switch (mat_id) {
+ case id_ANEOS_iron:
+ printf(" ANEOS iron \n");
+ break;
+
+ case id_MANEOS_forsterite:
+ printf(" MANEOS forsterite \n");
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d \n", mat_id);
+ };
+ break;
+
+ // SESAME
+ case type_SESAME:
+ switch (mat_id) {
+ case id_SESAME_iron:
+ printf(" SESAME iron \n");
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d \n", mat_id);
+ };
+ break;
+
+ default:
+ error("Unknown material type! mat_id = %d \n", mat_id);
+ }
+
+ // Convert to internal units (Earth masses and radii)
+ units_init(&us, 5.9724e27, 6.3710e8, 1.f, 1.f, 1.f);
+ log_rho_min -= logf(units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY));
+ log_rho_max -= logf(units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY));
+ log_u_min +=
+ logf(J_kg_to_erg_g /
+ units_cgs_conversion_factor(&us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
+ log_u_max +=
+ logf(J_kg_to_erg_g /
+ units_cgs_conversion_factor(&us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
+
+ // Initialise the EOS materials
+ eos_init(&eos, phys_const, &us, params);
+
+ // Output file
+ sprintf(filename, "testEOS_rho_u_P_%d.txt", mat_id);
+ FILE *f = fopen(filename, "w");
+ if (f == NULL) {
+ printf("Could not open output file!\n");
+ exit(EXIT_FAILURE);
+ }
+
+ if (do_output == 1) {
+ fprintf(f, "Density Sp.Int.Energy mat_id \n");
+ fprintf(f, "%d %d %d \n", num_rho, num_u, mat_id);
+ }
+
+ // Densities
+ log_rho = log_rho_min;
+ for (int i = 0; i < num_rho; i++) {
+ A1_rho[i] = exp(log_rho);
+ log_rho += log_rho_step;
+
+ if (do_output == 1)
+ fprintf(f, "%.6e ",
+ A1_rho[i] * units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY));
+ }
+ if (do_output == 1) fprintf(f, "\n");
+
+ // Sp. int. energies
+ log_u = log_u_min;
+ for (int i = 0; i < num_u; i++) {
+ A1_u[i] = exp(log_u);
+ log_u += log_u_step;
+
+ if (do_output == 1)
+ fprintf(f, "%.6e ", A1_u[i] *
+ units_cgs_conversion_factor(&us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
+ }
+ if (do_output == 1) fprintf(f, "\n");
+
+ // Pressures
+ for (int i = 0; i < num_rho; i++) {
+ rho = A1_rho[i];
+
+ for (int j = 0; j < num_u; j++) {
+ P = gas_pressure_from_internal_energy(rho, A1_u[j], mat_id);
+
+ if (do_output == 1)
+ fprintf(f, "%.6e ",
+ P * units_cgs_conversion_factor(&us, UNIT_CONV_PRESSURE));
+ }
+
+ if (do_output == 1) fprintf(f, "\n");
+ }
+ fclose(f);
+
+ return 0;
+}
+#else
+int main() {
+ return 0;
+}
+#endif
diff --git a/tests/testEOS.py b/tests/testEOS.py
new file mode 100644
index 0000000000000000000000000000000000000000..363bab200b58c65fa24cc033af4b8d3c04b7b503
--- /dev/null
+++ b/tests/testEOS.py
@@ -0,0 +1,176 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 2018 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+"""
+Plot the output of testEOS to show how the equation of state pressure varies
+with density and specific internal energy.
+
+Usage:
+ python testEOS.py (mat_id)
+
+ Sys args (optional):
+ mat_id | int | Material ID, see equation_of_state.h for the options.
+ Default: 201 (= HM80_ice).
+
+Text file contains:
+ header
+ num_rho num_u mat_id # Header info
+ rho_0 rho_1 rho_2 ... rho_num_rho # Array of densities, rho
+ u_0 u_1 u_2 ... u_num_u # Array of energies, u
+ P_0_0 P_0_1 ... P_0_num_u # Array of pressures, P(rho, u)
+ P_1_0 ... ... P_1_num_u
+ ... ... ... ...
+ P_num_rho_0 ... P_num_rho_num_u
+
+Note that the values tested extend beyond the range that most EOS are
+designed for (e.g. outside table limits), to help test the EOS in case of
+unexpected particle behaviour.
+"""
+
+# ========
+# Modules and constants
+# ========
+from __future__ import division
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import sys
+
+# Material types (copied from src/equation_of_state/planetary/equation_of_state.h)
+type_factor = 100
+Di_type = {
+ 'Til' : 1,
+ 'HM80' : 2,
+ 'ANEOS' : 3,
+ 'SESAME' : 4,
+}
+Di_material = {
+ # Tillotson
+ 'Til_iron' : Di_type['Til']*type_factor,
+ 'Til_granite' : Di_type['Til']*type_factor + 1,
+ 'Til_water' : Di_type['Til']*type_factor + 2,
+ # Hubbard & MacFarlane (1980) Uranus/Neptune
+ 'HM80_HHe' : Di_type['HM80']*type_factor, # Hydrogen-helium atmosphere
+ 'HM80_ice' : Di_type['HM80']*type_factor + 1, # H20-CH4-NH3 ice mix
+ 'HM80_rock' : Di_type['HM80']*type_factor + 2, # SiO2-MgO-FeS-FeO rock mix
+ # ANEOS
+ 'ANEOS_iron' : Di_type['ANEOS']*type_factor,
+ 'MANEOS_forsterite' : Di_type['ANEOS']*type_factor + 1,
+ # SESAME
+ 'SESAME_iron' : Di_type['SESAME']*type_factor,
+}
+# Invert so the mat_id are the keys
+Di_mat_id = {mat_id : mat for mat, mat_id in Di_material.iteritems()}
+
+# Unit conversion
+Ba_to_Mbar = 1e-12
+erg_g_to_J_kg = 1e-4
+
+if __name__ == '__main__':
+ # Sys args
+ try:
+ mat_id = int(sys.argv[1])
+ except IndexError:
+ mat_id = Di_material['HM80_ice']
+
+ # Check the material
+ try:
+ mat = Di_mat_id[mat_id]
+ print mat
+ sys.stdout.flush()
+ except KeyError:
+ print "Error: unknown material ID! mat_id = %d" % mat_id
+ print "Materials:"
+ for mat_id, mat in sorted(Di_mat_id.iteritems()):
+ print " %s%s%d" % (mat, (20 - len("%s" % mat))*" ", mat_id)
+
+ filename = "testEOS_rho_u_P_%d.txt" % mat_id
+
+ # Load the header info and density and energy arrays
+ with open(filename) as f:
+ f.readline()
+ num_rho, num_u, mat_id = np.array(f.readline().split(), dtype=int)
+ A1_rho = np.array(f.readline().split(), dtype=float)
+ A1_u = np.array(f.readline().split(), dtype=float)
+
+ # Load pressure array
+ A2_P = np.loadtxt(filename, skiprows=4)
+
+ # Convert pressures from cgs Barye to Mbar
+ A2_P *= Ba_to_Mbar
+ # Convert energies from cgs to SI
+ A1_u *= erg_g_to_J_kg
+
+ # Check that the numbers are right
+ assert A1_rho.shape == (num_rho,)
+ assert A1_u.shape == (num_u,)
+ assert A2_P.shape == (num_rho, num_u)
+
+ # Plot
+ plt.figure(figsize=(7, 7))
+ ax = plt.gca()
+
+ # P(rho) at fixed u
+ num_u_fix = 9
+ A1_idx = np.floor(np.linspace(0, num_u - 1, num_u_fix)).astype(int)
+ A1_colour = matplotlib.cm.rainbow(np.linspace(0, 1, num_u_fix))
+
+ for i, idx in enumerate(A1_idx):
+ plt.plot(A1_rho, A2_P[:, idx], c=A1_colour[i],
+ label=r"%.2e" % A1_u[idx])
+
+ plt.legend(title="Sp. Int. Energy (J kg$^{-1}$)")
+ plt.xscale('log')
+ plt.yscale('log')
+ plt.xlabel(r"Density (g cm$^{-3}$)")
+ plt.ylabel(r"Pressure (Mbar)")
+ plt.title(mat)
+ plt.tight_layout()
+
+ plt.savefig("testEOS_%d.png" % mat_id)
+ plt.close()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testEOS.sh b/tests/testEOS.sh
new file mode 100755
index 0000000000000000000000000000000000000000..411ac746be186bfe5758e03c2a852e081daefd10
--- /dev/null
+++ b/tests/testEOS.sh
@@ -0,0 +1,23 @@
+#!/bin/bash
+
+echo ""
+
+rm -f testEOS_rho_u_P_*.txt
+
+echo "Running testEOS for each planetary material"
+
+A1_mat_id=(
+ 100
+ 101
+ 102
+ 200
+ 201
+ 202
+)
+
+for mat_id in "${A1_mat_id[@]}"
+do
+ ./testEOS "$mat_id" 1
+done
+
+exit $?
diff --git a/tests/testEOS_plot.sh b/tests/testEOS_plot.sh
new file mode 100755
index 0000000000000000000000000000000000000000..39108c5e19d8f4474de508e205951a1fd0aebcc9
--- /dev/null
+++ b/tests/testEOS_plot.sh
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+echo ""
+
+echo "Plotting testEOS output for each planetary material"
+
+A1_mat_id=(
+ 100
+ 101
+ 102
+ 200
+ 201
+ 202
+)
+
+for mat_id in "${A1_mat_id[@]}"
+do
+ python ./testEOS.py "$mat_id"
+done
+
+exit $?