# Source code for burnman.eos.murnaghan

# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
# the Earth and Planetary Sciences
# Copyright (C) 2012 - 2021 by the BurnMan team, released under the GNU
# GPL v2 or later.

from . import equation_of_state as eos
import warnings
import numpy as np

def volume(pressure, V_0, K_0, Kprime_0):
return V_0 * np.power(1. + (pressure * Kprime_0 / K_0), -1./Kprime_0)

def pressure(volume, V_0, K_0, Kprime_0):
return K_0 / Kprime_0 * (np.power(volume / V_0, -Kprime_0) - 1.)

def bulk_modulus(pressure, K_0, Kprime_0):
return K_0 + pressure * Kprime_0

def energy(volume, E_0, V_0, K_0, Kprime_0):
Vrel = volume / V_0
return (E_0 + K_0*V_0 * (np.power(Vrel, 1. - Kprime_0)
/ (Kprime_0*(Kprime_0 - 1))
+ Vrel / Kprime_0 - 1./(Kprime_0 - 1.)))

def intVdP(pressure, V_0, K_0, Kprime_0):
return (V_0 * K_0
* ((np.power(1. + (pressure * Kprime_0 / K_0),
1.-(1./Kprime_0))) - 1.)
/ (Kprime_0 - 1.))

[docs]class Murnaghan(eos.EquationOfState):

"""
Base class for the isothermal Murnaghan equation of state,
as described in :cite:Murnaghan1944.
"""

[docs]    def volume(self, pressure, temperature, params):
"""
Returns volume :math:[m^3] as a function of pressure :math:[Pa].
"""
return volume(pressure,
params['V_0'], params['K_0'], params['Kprime_0'])

[docs]    def pressure(self, temperature, volume, params):
return pressure(volume,
params['V_0'], params['K_0'], params['Kprime_0'])

[docs]    def isothermal_bulk_modulus(self, pressure, temperature, volume, params):
"""
Returns isothermal bulk modulus :math:K_T :math:[Pa]
as a function of pressure :math:[Pa],
temperature :math:[K] and volume :math:[m^3].
"""
return bulk_modulus(pressure, params['K_0'], params['Kprime_0'])

[docs]    def adiabatic_bulk_modulus(self, pressure, temperature, volume, params):
"""
Returns adiabatic bulk modulus :math:K_s of the mineral. :math:[Pa].
"""
return bulk_modulus(pressure, params['K_0'], params['Kprime_0'])

[docs]    def shear_modulus(self, pressure, temperature, volume, params):
"""
Returns shear modulus :math:G of the mineral. :math:[Pa]
Currently not included in the Murnghan EOS, so omitted.
"""
return 0.

[docs]    def entropy(self, pressure, temperature, volume, params):
"""
Returns the molar entropy :math:\mathcal{S} of the mineral.
:math:[J/K/mol]
"""
return 0.

[docs]    def molar_internal_energy(self, pressure, temperature, volume, params):
"""
Returns the internal energy :math:\mathcal{E} of the mineral.
:math:[J/mol]
"""
return energy(volume,
params['E_0'], params['V_0'],
params['K_0'], params['Kprime_0'])

[docs]    def gibbs_free_energy(self, pressure, temperature, volume, params):
"""
Returns the Gibbs free energy :math:\mathcal{G} of the mineral.
:math:[J/mol]
"""
# G = E + PV
return (self.molar_internal_energy(pressure, temperature,
volume, params)
+ volume*pressure)

[docs]    def molar_heat_capacity_v(self, pressure, temperature, volume, params):
"""
Since this equation of state does not contain temperature effects,
return a very large number. :math:[J/K/mol]
"""
return 1.e99

[docs]    def molar_heat_capacity_p(self, pressure, temperature, volume, params):
"""
Since this equation of state does not contain temperature effects,
return a very large number. :math:[J/K/mol]
"""
return 1.e99

[docs]    def thermal_expansivity(self, pressure, temperature, volume, params):
"""
Since this equation of state does not contain temperature effects,
return zero. :math:[1/K]
"""
return 0.

[docs]    def grueneisen_parameter(self, pressure, temperature, volume, params):
"""
Since this equation of state does not contain temperature effects,
return zero. :math:[unitless]
"""
return 0.

[docs]    def validate_parameters(self, params):
"""
Check for existence and validity of the parameters
"""

if 'E_0' not in params:
params['E_0'] = 0.
if 'P_0' not in params:
params['P_0'] = 0.

# G is not included in the Murnaghan EOS so we shall set them to NaN's
if 'G_0' not in params:
params['G_0'] = float('nan')
if 'Gprime_0' not in params:
params['Gprime_0'] = float('nan')

# check that all the required keys are in the dictionary
expected_keys = ['V_0', 'K_0', 'Kprime_0']
for k in expected_keys:
if k not in params:
raise KeyError('params object missing parameter : ' + k)

# now check that the values are reasonable.  I mostly just
# made up these values from experience, and we are only
# raising a warning.  Better way to do this? [IR]
if params['V_0'] < 1.e-7 or params['V_0'] > 1.e-3:
warnings.warn('Unusual value for V_0', stacklevel=2)
if params['K_0'] < 1.e9 or params['K_0'] > 1.e13:
warnings.warn('Unusual value for K_0', stacklevel=2)
if params['Kprime_0'] < -5. or params['Kprime_0'] > 30.:
warnings.warn('Unusual value for Kprime_0', stacklevel=2)