testEOS.py 6.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
###############################################################################
 # 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/>.
 #
 ##############################################################################
"""
21 22
Plot the output of testEOS to show how the equation of state pressure and sound
speed varies with density and specific internal energy.
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39

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
40 41 42 43
    c_0_0   c_0_1   ...     c_0_num_u           # Array of sound speeds, c(rho, u)
    c_1_0   ...     ...     c_1_num_u
    ...     ...     ...     ...
    c_num_rho_0     ...     c_num_rho_num_u
44 45

Note that the values tested extend beyond the range that most EOS are
46
designed for (e.g. outside table limits), to help test the EOS in cases of
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
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,
65
    'SESAME'    : 3,
66 67 68 69 70 71 72 73 74 75 76
}
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
    # SESAME
77 78 79 80
    'SESAME_iron'   : Di_type['SESAME']*type_factor,        # 2140
    'SESAME_basalt' : Di_type['SESAME']*type_factor + 1,    # 7530
    'SESAME_water'  : Di_type['SESAME']*type_factor + 2,    # 7154
    'SS08_water'    : Di_type['SESAME']*type_factor + 3,    # Senft & Stewart (2008)
81 82 83 84 85 86 87
}
# 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
88
cm_s_to_m_s = 1e-2
89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107

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)

108
    filename = "testEOS_rho_u_P_c_%d.txt" % mat_id
109 110 111 112 113 114 115 116

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

117
    # Load pressure and soundspeed arrays
118
    A2_P = np.loadtxt(filename, skiprows=4)
119 120
    A2_c = A2_P[num_rho:]
    A2_P = A2_P[:num_rho]
121 122 123

    # Convert energies from cgs to SI
    A1_u *= erg_g_to_J_kg
124 125 126 127
    # Convert pressures from cgs (Barye) to Mbar
    A2_P *= Ba_to_Mbar
    # Convert sound speeds from cgs to SI
    A1_u *= cm_s_to_m_s
128 129 130 131 132

    # 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)
133
    assert A2_c.shape == (num_rho, num_u)
134 135

    # Plot
136
    # Pressure: P(rho) at fixed u
137 138 139
    plt.figure(figsize=(7, 7))
    ax = plt.gca()

140
    A1_colour = matplotlib.cm.rainbow(np.linspace(0, 1, num_u))
141

142 143 144 145 146 147
    for i_u, u in enumerate(A1_u):
        if i_u%10 == 0:
            plt.plot(A1_rho, A2_P[:, i_u], c=A1_colour[i_u],
                     label=r"%.2e" % u)
        else:
            plt.plot(A1_rho, A2_P[:, i_u], c=A1_colour[i_u])
148 149 150 151 152 153 154 155 156

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

157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181
    plt.savefig("testEOS_P_%d.png" % mat_id)
    plt.close()

    # Sound speed: c(rho) at fixed u
    plt.figure(figsize=(7, 7))
    ax = plt.gca()

    A1_colour = matplotlib.cm.rainbow(np.linspace(0, 1, num_u))

    for i_u, u in enumerate(A1_u):
        if i_u%10 == 0:
            plt.plot(A1_rho, A2_c[:, i_u], c=A1_colour[i_u],
                     label=r"%.2e" % u)
        else:
            plt.plot(A1_rho, A2_c[:, i_u], c=A1_colour[i_u])

    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"Sound Speed (m s^{-1})")
    plt.title(mat)
    plt.tight_layout()

    plt.savefig("testEOS_c_%d.png" % mat_id)
182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
    plt.close()