Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix BM4, MT pressure calcs #611

Merged
merged 2 commits into from
Nov 25, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions burnman/calibrants/Decker_1971.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@

from burnman.eos.mie_grueneisen_debye import MGDBase
from burnman.classes.calibrant import Calibrant
from burnman.eos.birch_murnaghan_4th import birch_murnaghan_fourth
from burnman.eos.birch_murnaghan_4th import pressure_birch_murnaghan_fourth


class NaCl_B1(Calibrant):
Expand All @@ -27,7 +27,7 @@ class NaCl_B1(Calibrant):
def __init__(self):
def _pressure_Decker_NaCl(volume, temperature, params):
# Isothermal pressure (GPa)
P0 = birch_murnaghan_fourth(params["V_0"] / volume, params)
P0 = pressure_birch_murnaghan_fourth(params["V_0"] / volume, params)

# Thermal pressure
thermal_model = MGDBase()
Expand Down
44 changes: 33 additions & 11 deletions burnman/eos/birch_murnaghan_4th.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,8 @@ def bulk_modulus_fourth(volume, params):
modulus. Pressure must be in :math:`[Pa]`.
"""

x = params["V_0"] / volume
f = 0.5 * (pow(x, 2.0 / 3.0) - 1.0)
invVrel = params["V_0"] / volume
f = 0.5 * (pow(invVrel, 2.0 / 3.0) - 1.0)

Xi = (3.0 / 4.0) * (4.0 - params["Kprime_0"])
Zeta = (3.0 / 8.0) * (
Expand All @@ -45,24 +45,46 @@ def bulk_modulus_fourth(volume, params):


def volume_fourth_order(pressure, params):
func = lambda x: birch_murnaghan_fourth(params["V_0"] / x, params) - pressure
"""
Volume according to the fourth-order
Birch-Murnaghan equation of state.

:param pressure: Pressure in the same units that are supplied for the reference bulk
modulus (params['K_0'])
:type pressure: float
:param params: parameter dictionary
:type params: dictionary
:return: V/V_0
:rtype: float
"""

def delta_pressure(x):
return pressure_birch_murnaghan_fourth(params["V_0"] / x, params) - pressure

try:
sol = bracket(func, params["V_0"], 1.0e-2 * params["V_0"])
except:
sol = bracket(delta_pressure, params["V_0"], 1.0e-2 * params["V_0"])
except ValueError:
raise ValueError(
"Cannot find a volume, perhaps you are outside of the range of validity for the equation of state?"
)
return opt.brentq(func, sol[0], sol[1])
return opt.brentq(delta_pressure, sol[0], sol[1])


def birch_murnaghan_fourth(x, params):
def pressure_birch_murnaghan_fourth(invVrel, params):
"""
equation for the fourth order birch-murnaghan equation of state, returns
pressure in the same units that are supplied for the reference bulk
Pressure according to the fourth-order
Birch-Murnaghan equation of state.

:param invVrel: V/V_0
:type invVrel: float or numpy array
:param params: parameter dictionary
:type params: dictionary
:return: Pressure in the same units that are supplied for the reference bulk
modulus (params['K_0'])
:rtype: float or numpy array
"""

f = 0.5 * (pow(x, 2.0 / 3.0) - 1.0)
f = 0.5 * (pow(invVrel, 2.0 / 3.0) - 1.0)
Xi = (3.0 / 4.0) * (4.0 - params["Kprime_0"])
Zeta = (3.0 / 8.0) * (
(params["K_0"] * params["Kprime_prime_0"])
Expand Down Expand Up @@ -93,7 +115,7 @@ def volume(self, pressure, temperature, params):
return volume_fourth_order(pressure, params)

def pressure(self, temperature, volume, params):
return birch_murnaghan_fourth(volume / params["V_0"], params)
return pressure_birch_murnaghan_fourth(params["V_0"] / volume, params)

def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params):
"""
Expand Down
22 changes: 14 additions & 8 deletions burnman/eos/modified_tait.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,15 +30,21 @@ def tait_constants(params):
return a, b, c


def modified_tait(x, params):
def pressure_modified_tait(Vrel, params):
"""
equation for the modified Tait equation of state, returns
pressure in the same units that are supplied for the reference bulk
modulus (params['K_0'])
Pressure according to the modified Tait equation of state.
EQ 2 from Holland and Powell, 2011

:param Vrel: V/V_0
:type Vrel: float or numpy array
:param params: Parameter dictionary
:type params: dictionary
:return: pressure in the same units that are supplied for the reference bulk
modulus (params['K_0'])
:rtype: float or numpy array
"""
a, b, c = tait_constants(params)
return (np.power((x + a - 1.0) / a, -1.0 / c) - 1.0) / b + params["P_0"]
return (np.power((Vrel + a - 1.0) / a, -1.0 / c) - 1.0) / b + params["P_0"]


def volume(pressure, params):
Expand All @@ -47,8 +53,8 @@ def volume(pressure, params):
EQ 12
"""
a, b, c = tait_constants(params)
x = 1 - a * (1.0 - np.power((1.0 + b * (pressure - params["P_0"])), -1.0 * c))
return x * params["V_0"]
Vrel = 1.0 - a * (1.0 - np.power((1.0 + b * (pressure - params["P_0"])), -1.0 * c))
return Vrel * params["V_0"]


def bulk_modulus(pressure, params):
Expand Down Expand Up @@ -113,7 +119,7 @@ def pressure(self, temperature, volume, params):
"""
Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
"""
return modified_tait(params["V_0"] / volume, params)
return pressure_modified_tait(volume / params["V_0"], params)

def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params):
"""
Expand Down
Loading