diff --git a/.gitignore b/.gitignore
index 1a2e83a9d8c62f48cfa9db9f45b1b40bec80b449..27fa7cf814b81dc5ba4bf697a8a4e0e6a1393333 100644
--- a/.gitignore
+++ b/.gitignore
@@ -111,6 +111,8 @@ tests/benchmarkInteractions
tests/testGravityDerivatives
tests/testPotentialSelf
tests/testPotentialPair
+tests/testEOS
+tests/testEOS*.txt
theory/latex/swift.pdf
theory/SPH/Kernels/kernels.pdf
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 643d89aaf05e31b6ea1c304c4fac6c28f7e89d03..cead0cc9cd16ae06a08e051c913db36d8cb43558 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -1,16 +1,16 @@
# This file is part of SWIFT.
# Copyright (c) 2015 matthieu.schaller@durham.ac.uk.
-#
+#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
-#
+#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
-#
+#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
@@ -27,7 +27,7 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
testMatrixInversion testThreadpool testDump testLogger testInteractions.sh \
testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
- testPotentialPair
+ testPotentialPair testEOS
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
@@ -37,7 +37,7 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testAdiabaticIndex testRiemannExact testRiemannTRRS \
testRiemannHLLC testMatrixInversion testDump testLogger \
testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
- testGravityDerivatives testPotentialSelf testPotentialPair
+ testGravityDerivatives testPotentialSelf testPotentialPair testEOS
# Rebuild tests when SWIFT is updated.
$(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -105,6 +105,8 @@ testPotentialSelf_SOURCES = testPotentialSelf.c
testPotentialPair_SOURCES = testPotentialPair.c
+testEOS_SOURCES = testEOS.c
+
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testActivePair.sh \
test27cells.sh test27cellsPerturbed.sh testParser.sh testPeriodicBC.sh \
diff --git a/tests/testEOS.c b/tests/testEOS.c
new file mode 100644
index 0000000000000000000000000000000000000000..aa04a819825ecf3ec257731ab6fec7b60f2629d4
--- /dev/null
+++ b/tests/testEOS.c
@@ -0,0 +1,241 @@
+/*******************************************************************************
+ * 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 "swift.h"
+#include "equation_of_state.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.
+ *
+ * So far only has the Hubbard & MacFarlane (1980) equations of state.
+ *
+ * Sys args:
+ * mat_id | int | Material ID, see equation_of_state.h for the options.
+ * Default: 21 (= 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 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
+ *
+ */
+
+int main(int argc, char *argv[]) {
+ int mat_id, do_output;
+ struct HM80_params mat;
+ float rho, log_rho, log_u, P;
+ int num_rho, num_u;
+ struct unit_system us;
+ const struct phys_const *phys_const = 0; // Unused placeholder
+ const struct swift_params *params = 0; // Unused placeholder
+
+ // Check the number of system arguments and set defaults if not provided
+ switch(argc) {
+ case 1:
+ // Default both
+ mat_id = HM80_ice;
+ do_output = 0;
+ break;
+
+ case 2:
+ // Read mat_id, default do_output
+ mat_id = atoi(argv[1]);
+ do_output = 0;
+ 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 = HM80_ice; // Ignored, just here to keep the compiler happy
+ do_output = 0;
+ };
+
+ /* Greeting message */
+ printf("This is %s\n", package_description());
+
+ // Select the material parameters
+ switch(mat_id) {
+ case HM80_HHe:
+ printf("HM80_HHe \n");
+ set_HM80_HHe(&mat);
+ load_HM80_table(&mat, HM80_HHe_table_file);
+ break;
+
+ case HM80_ice:
+ printf("HM80_ice \n");
+ set_HM80_ice(&mat);
+ load_HM80_table(&mat, HM80_ice_table_file);
+ break;
+
+ case HM80_rock:
+ printf("HM80_rock \n");
+ set_HM80_rock(&mat);
+ load_HM80_table(&mat, HM80_rock_table_file);
+ break;
+
+ default:
+ error("Unknown material ID! mat_id = %d", mat_id);
+ set_HM80_rock(&mat); // Ignored, just here to keep the compiler happy
+ load_HM80_table(&mat, HM80_rock_table_file);
+ };
+
+ // Convert to internal units
+ units_init(&us, 5.9724e27, 6.3710e8, 1, 1, 1);
+ convert_units_HM80(&mat, &us);
+
+ eos_init(&eos, phys_const, &us, params);
+
+ // Output file
+ FILE *f = fopen("testEOS_rho_u_P.txt", "w");
+ if (f == NULL) {
+ printf("Could not open output file!\n");
+ exit(EXIT_FAILURE);
+ }
+
+ num_rho = (mat.log_rho_max - mat.log_rho_min) / mat.log_rho_step;
+ num_u = (mat.log_u_max - mat.log_u_min) / mat.log_u_step;
+ 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);
+ }
+
+ // Arrays of densities and energies
+ float A1_rho[num_rho], A1_u[num_u];
+
+ log_rho = mat.log_rho_min;
+ for (int i = 0; i < num_rho; i++) {
+ A1_rho[i] = exp(log_rho);
+ log_rho += mat.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");
+ log_u = mat.log_u_min;
+ for (int i = 0; i < num_u; i++) {
+ A1_u[i] = exp(log_u);
+ log_u += mat.log_u_step;
+
+ if (do_output == 1)
+ fprintf(f, "%.6e ", A1_u[i] *
+ units_cgs_conversion_factor(&us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
+ }
+
+ // Pressures P(rho, u)
+ if (do_output == 1) fprintf(f, "\n");
+ 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.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;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/testEOS.sh b/tests/testEOS.sh
new file mode 100644
index 0000000000000000000000000000000000000000..4ffb43eecf03770f87220f4ee5a0d18e52ea5cf3
--- /dev/null
+++ b/tests/testEOS.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+
+echo ""
+
+rm -f testEOS_rho_u_P.txt
+
+echo "Running ./testEOS 21 1 0"
+
+./testEOS 21 1 0
+
+exit $?
diff --git a/tests/testEOS_plot.py b/tests/testEOS_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..83b512683ed9c2a0027c8792a0ef3c4c269b365e
--- /dev/null
+++ b/tests/testEOS_plot.py
@@ -0,0 +1,153 @@
+###############################################################################
+ # 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.
+
+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
+"""
+
+# ========
+# Modules and constants
+# ========
+from __future__ import division
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+
+filename = "testEOS_rho_u_P.txt"
+
+# Material types (copied from src/equation_of_state/planetary/equation_of_state.h)
+type_factor = 10
+Di_type = {
+ 'Till' : 1,
+ 'HM80' : 2,
+ 'ANEOS' : 3,
+ 'SESAME' : 4,
+}
+Di_material = {
+ # Tillotson
+ 'Til_iron' : Di_type['Till']*type_factor,
+ 'Til_granite' : Di_type['Till']*type_factor + 1,
+ 'Til_water' : Di_type['Till']*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()}
+
+Ba_to_Mbar = 1e-12
+erg_g_to_J_kg = 1e-4
+
+# ========
+# Main
+# ========
+# 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)
+try:
+ mat = Di_mat_id[mat_id]
+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)
+
+# 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.show()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+