assert Phi=0 at standard state

geodynamics · Mar 5, 2024 · c9746ba · c9746ba
1 parent f25df4a
commit c9746ba
Show file tree

Hide file tree

Showing 5 changed files with 18 additions and 12 deletions.
diff --git a/burnman/classes/anisotropicmineral.py b/burnman/classes/anisotropicmineral.py
@@ -204,6 +204,12 @@ def __init__(
  self.anisotropic_params = anisotropic_parameters
  self.psi_function = psi_function
 
+ Psi_Voigt0 = self.psi_function(0., 0., self.anisotropic_params)[0]
+ if not np.all(Psi_Voigt0 == 0.):
+ raise ValueError("All elements of Psi should evaluate to zero at "
+ "standard state. The current array evaluates to: "
+ f"{Psi_Voigt0}")
+
  # cell_vectors is the transpose of the cell tensor M
  self.cell_vectors_0 = cell_parameters_to_vectors(cell_parameters)
 

diff --git a/contrib/anisotropic_eos/cubic_fitting_second_form.py b/contrib/anisotropic_eos/cubic_fitting_second_form.py
@@ -61,8 +61,8 @@ def psi_func(f, Pth, params):
  Psi = (
  0.0
  + params["a"] * f
- + params["b_1"] * np.exp(params["c_1"] * f)
- + params["b_2"] * np.exp(params["c_2"] * f)
+ + params["b_1"] * (np.exp(params["c_1"] * f) - 1.)
+ + params["b_2"] * (np.exp(params["c_2"] * f) - 1.)
  )
  Psi = Psi + Pth / 1.0e9 * (
  params["b_3"] * np.exp(params["c_3"] * f)
@@ -379,8 +379,8 @@ def cubic_misfit(x, imin):
  f'the isotropic model: $V_0$: {m.params["V_0"]*1.e6:.5f} cm$^3$/mol, '
  f'$K_0$: {m.params["K_0"]/1.e9:.5f} GPa, '
  f'$K\'_0$: {m.params["Kprime_0"]:.5f}, '
- f'$\Theta_0$: {m.params["Debye_0"]:.5f} K, '
- f'$\gamma_0$: {m.params["grueneisen_0"]:.5f}, '
+ f'$\\Theta_0$: {m.params["Debye_0"]:.5f} K, '
+ f'$\\gamma_0$: {m.params["grueneisen_0"]:.5f}, '
  f'and $q_0$: {m.params["q_0"]:.5f}.'
 )
 

diff --git a/contrib/anisotropic_eos/cubic_fitting_second_form_atherm.py b/contrib/anisotropic_eos/cubic_fitting_second_form_atherm.py
@@ -53,8 +53,8 @@ def psi_func(f, Pth, params):
  Psi = (
  0.0
  + params["a"] * f
- + params["b_1"] * np.exp(params["c_1"] * f)
- + params["b_2"] * np.exp(params["c_2"] * f)
+ + params["b_1"] * (np.exp(params["c_1"] * f) - 1.)
+ + params["b_2"] * (np.exp(params["c_2"] * f) - 1.)
  )
  dPsidPth = 0.0 * params["b_1"]
  return (Psi, dPsidf, dPsidPth)

diff --git a/contrib/anisotropic_eos/cubic_fitting_second_form_only_therm.py b/contrib/anisotropic_eos/cubic_fitting_second_form_only_therm.py
@@ -114,8 +114,8 @@ def psi_func(f, Pth, params):
  Psi = (
  0.0
  + params["a"] * f
- + params["b_1"] * np.exp(params["c_1"] * f)
- + params["b_2"] * np.exp(params["c_2"] * f)
+ + params["b_1"] * (np.exp(params["c_1"] * f) - 1.)
+ + params["b_2"] * (np.exp(params["c_2"] * f) - 1.)
  )
  Psi = Psi + Pth / 1.0e9 * (
  params["b_3"] * np.exp(params["c_3"] * f)

diff --git a/contrib/anisotropic_eos/orthorhombic_fitting_second_form.py b/contrib/anisotropic_eos/orthorhombic_fitting_second_form.py
@@ -100,8 +100,8 @@ def psi_func(f, Pth, params):
  Psi = (
  0.0
  + params["a"] * f
- + params["b_1"] * np.exp(params["c_1"] * f)
- + params["b_2"] * np.exp(params["c_2"] * f)
+ + params["b_1"] * (np.exp(params["c_1"] * f) - 1.)
+ + params["b_2"] * (np.exp(params["c_2"] * f) - 1.)
  + params["d"] * Pth / 1.0e9
  )
  dPsidPth = params["d"] / 1.0e9
@@ -393,8 +393,8 @@ def orthorhombic_misfit(x, imin):
  f'the isotropic model: $V_0$: {m.params["V_0"]*1.e6:.5f} cm$^3$/mol, '
  f'$K_0$: {m.params["K_0"]/1.e9:.5f} GPa, '
  f'$K\'_0$: {m.params["Kprime_0"]:.5f}, '
- f'$\Theta_0$: {m.params["Debye_0"]:.5f} K, '
- f'$\gamma_0$: {m.params["grueneisen_0"]:.5f}, '
+ f'$\\Theta_0$: {m.params["Debye_0"]:.5f} K, '
+ f'$\\gamma_0$: {m.params["grueneisen_0"]:.5f}, '
  f'and $q_0$: {m.params["q_0"]:.5f}.'
  )