diff --git a/tests/testEOS.c b/tests/testEOS.c deleted file mode 100644 index 097e4db451455d4a4e72ebdd72faf0ee729fe5bc..0000000000000000000000000000000000000000 --- a/tests/testEOS.c +++ /dev/null @@ -1,201 +0,0 @@ -/******************************************************************************* - * 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. - * - * 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 deleted file mode 100644 index 4ffb43eecf03770f87220f4ee5a0d18e52ea5cf3..0000000000000000000000000000000000000000 --- a/tests/testEOS.sh +++ /dev/null @@ -1,11 +0,0 @@ -#!/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 deleted file mode 100644 index 83b512683ed9c2a0027c8792a0ef3c4c269b365e..0000000000000000000000000000000000000000 --- a/tests/testEOS_plot.py +++ /dev/null @@ -1,153 +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/>. - # - ############################################################################## -""" -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() - - - - - - - - - - - - - - - - - - - - - - - - - - - -