Skip to content
Snippets Groups Projects
Commit b3f49060 authored by Jacob Kegerreis's avatar Jacob Kegerreis
Browse files

Remove testEOS for now

parent 8f6cd1ca
Branches
Tags
1 merge request!545Add support for equations of state related to planetary physics
Hide whitespace changes
Inline Side-by-side
tests/testEOS.c deleted 100644 → 0
+ 0
− 201
View file @ 8f6cd1ca
/*******************************************************************************
* 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;
}
tests/testEOS.sh deleted 100644 → 0
+ 0
− 11
View file @ 8f6cd1ca
#!/bin/bash
echo ""
rm -f testEOS_rho_u_P.txt
echo "Running ./testEOS 21 1 0"
./testEOS 21 1 0
exit $?
tests/testEOS_plot.py deleted 100644 → 0
+ 0
− 153
View file @ 8f6cd1ca
###############################################################################
# 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()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment