Add back improved testEOS

.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

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

tests/testEOS.c

+/*******************************************************************************
+ * 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

tests/testEOS.py

+###############################################################################
+ # 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()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

tests/testEOS.sh

+#!/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 $?

tests/testEOS_plot.sh

+#!/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 $?