Skip to content

Commit

Permalink
Merge pull request #588 from bobmyhill/anisotropic_feldspar
Browse files Browse the repository at this point in the history 
anisotropic plagioclase
  • Loading branch information
@bobmyhill
bobmyhill committed Mar 27, 2024
2 parents c2f447f + acb682e commit dd46dcb
Show file tree
Hide file tree
Showing 9 changed files with 1,135 additions and 0 deletions.
34 changes: 34 additions & 0 deletions contrib/anisotropic_plagioclase/Readme.md
View file
@@ -0,0 +1,34 @@
## Companion code to
# An anisotropic equation of state for solid solutions
## R. Myhill

The python scripts contained in this directory accompany the paper
submitted by Myhill (2024).
There are comments in each of the scripts, so it is hoped that by reading
these in tandem with the original paper, the reader will get a feel
for how to use the code for their own purposes.

The scripts contained in this directory are as follows:

plagioclase_data.py
-------------------
This script packages all of the data in the data directory into convenient
dictionaries.

plagioclase_parameters.py
-------------------------
This script provides the final fitted parameters for the model,
along with the bounds used during fitting.

plagioclase_model.py
--------------------
This file makes the burnman solution objects from the parameters.

plagioclase_fit_eos.py
----------------------
This script allows the user to refine the model parameters
based on the experimental data.

plagioclase_model_plots.py
--------------------------
This script creates all of the plots presented in the paper.
14 changes: 14 additions & 0 deletions contrib/anisotropic_plagioclase/data/plagioclase_cell_tensor_Brown_2016.dat
View file
@@ -0,0 +1,14 @@
Sample p_an Space_group Z a a_err b b_err c c_err alpha alpha_err beta beta_err gamma gamma_err V V_err rho rho_err
An0 0.00 C1 4 8.1366 0.0002 12.7857 0.0002 7.1582 0.0003 94.253 0.002 116.605 0.002 87.756 0.002 663.98 0.03 2.623 0.003
An25 0.25 C1 4 8.1605 0.0015 12.8391 0.0006 7.1288 0.0003 93.833 0.013 116.440 0.005 89.124 0.005 667.20 0.13 2.653 0.003
An37 0.37 C1 4 8.16577 0.00009 12.8562 0.0001 7.11418 0.00009 93.622 0.001 116.278 0.001 89.679 0.001 668.123 0.012 2.666 0.003
An48 0.48 C1 4 8.1744 0.0002 12.8638 0.0003 7.1102 0.0002 93.525 0.003 116.236 0.001 89.915 0.003 669.10 0.03 2.683 0.003
An60 0.60 I1 8 8.1717 0.0002 12.8752 0.0002 14.2074 0.0003 93.4530 0.0011 116.078 0.001 91.4250 0.0011 1337.98 0.05 2.702 0.003
An67 0.67 I1 8 8.173 0.001 12.869 0.001 15.200 0.001 93.46 0.01 116.06 0.01 90.54 0.01 1338.1 0.1 2.721 0.003
An78 0.78 I1 8 8.1798 0.0003 12.8761 0.0003 14.1974 0.0006 93.423 0.003 116.054 0.003 90.705 0.003 1339.74 0.08 2.725 0.003
An96 0.96 P1 8 8.1789 0.0003 12.8717 0.0006 14.1765 0.0007 93.194 0.005 115.893 0.003 91.195 0.003 1338.84 0.1 2.757 0.003
# Table 1.Sample Chemistry, Unit Cell Parameters, and Densities Determined From Volume and Chemical Analyses
# Data for An0 are from Brown et al.[2006].
# Uncertainties in the last digits are given in parentheses.
# Journal of Geophysical Research: Solid Earth10.1002/2015JB012736
# BROWN ET AL.ELASTICITY OF PLAGIOCLASE FELDSPARS
11 changes: 11 additions & 0 deletions contrib/anisotropic_plagioclase/data/plagioclase_compressibilities_Brown_2016.dat
View file
@@ -0,0 +1,11 @@
# p_an beta_1 unc_beta_1 beta_2 unc_beta_2 beta_3 unc_beta_3 beta_4 unc_beta_4 beta_5 unc_beta_5 beta_6 unc_beta_6
0 0.00000e+00 1.10980e-02 1.63000e-04 3.37925e-03 9.98500e-05 3.65515e-03 1.35150e-04 2.72725e-04 1.01605e-04 9.43280e-04 1.00120e-04 2.05795e-03 9.98500e-05
1 2.00000e+01 8.69025e-03 3.80500e-05 3.28985e-03 1.00250e-04 3.89970e-03 9.98000e-05 -3.42245e-04 9.62550e-05 -1.22356e-04 1.02334e-04 1.94940e-03 9.92000e-05
2 3.70000e+01 7.32605e-03 8.69500e-05 3.18960e-03 1.00600e-04 3.95025e-03 9.98500e-05 -8.52945e-04 1.09625e-04 -5.16130e-04 1.00110e-04 1.75035e-03 9.98500e-05
3 4.60000e+01 6.54350e-03 1.08700e-04 3.12965e-03 1.00250e-04 3.97990e-03 9.90000e-05 -1.29145e-03 9.89500e-05 -4.22690e-04 1.00110e-04 1.42210e-03 1.00500e-04
4 4.80000e+01 7.00000e-03 1.19600e-04 3.19120e-03 5.11000e-05 3.88685e-03 1.13450e-04 -8.71655e-04 1.01605e-04 -3.58175e-04 1.02335e-04 1.71040e-03 9.99000e-05
5 6.00000e+01 6.17530e-03 1.09100e-04 3.21050e-03 8.05000e-05 4.02040e-03 1.00200e-04 -1.87435e-03 9.89500e-05 -9.45480e-04 1.02320e-04 2.11915e-03 9.98500e-05
6 6.80000e+01 5.71080e-03 1.06000e-04 3.27980e-03 3.06000e-05 4.06050e-03 1.01000e-04 -2.31285e-03 1.04250e-04 -1.26585e-03 1.02350e-04 2.26830e-03 9.99000e-05
7 7.80000e+01 5.43025e-03 5.30759e-05 3.27903e-03 5.12428e-05 4.05908e-03 5.94481e-05 -2.75935e-03 8.31911e-05 -1.51612e-03 9.74985e-05 2.61682e-03 8.07369e-05
8 8.90000e+01 5.03125e-03 9.97500e-05 3.34015e-03 4.99500e-05 4.05125e-03 1.00650e-04 -3.55345e-03 9.89500e-05 -2.09790e-03 9.79000e-05 2.28695e-03 9.98500e-05
9 9.60000e+01 4.91270e-03 1.43400e-04 3.06925e-03 1.00250e-04 3.93020e-03 9.98000e-05 -3.50000e-03 9.89000e-05 -1.29920e-03 1.98000e-04 2.14780e-03 9.92000e-05
26 changes: 26 additions & 0 deletions contrib/anisotropic_plagioclase/data/plagioclase_stiffness_tensor_Brown_2016.dat
View file
@@ -0,0 +1,26 @@
row M_an0 M_unc_an0 M_an25 M_unc_an25 M_an37 M_unc_an37 M_an48 M_unc_an48 M_an60 M_unc_an60 M_an67 M_unc_an67 M_an78 M_unc_an78 M_an96 M_unc_an96
values M M_unc M M_unc M M_unc M M_unc M M_unc M M_unc M M_unc M M_unc
p_an 0 0 25 25 37 37 48 48 60 60 67 67 78 78 96 96
# “Orthorhombic” Moduli
C11 68.3 0.8 87.1 1.3 96.2 1.6 104.6 1.9 109.3 1.7 120.3 4.2 120.4 2.6 132.2 3
C22 184.3 4.9 174.9 5.2 189.4 4.9 201.4 6.6 185.5 2.3 193.5 4.4 191.6 6.3 200.2 5.4
C33 180 3 166.1 4.7 171.9 4.5 172.8 5.1 164.1 1.9 171.9 5 163.7 5 163.9 4.1
C44 25 0.1 22.9 0.2 23.6 0.1 22.9 0.1 22.2 0.1 24 0.1 23.3 0.1 24.6 0.1
C55 26.9 0.1 29 0.2 33.1 0.3 33 0.3 33.1 0.2 35.5 0.3 32.8 0.3 36.6 0.2
C66 33.6 0.2 35 0.3 35.5 0.3 35.6 0.2 36.8 0.3 37.3 0.4 35 0.5 36 0.3
C12 32.2 1.6 43.9 2 46.1 2.5 51.5 2.8 53.1 1.1 54.4 3.7 56.6 3.4 64 3.5
C13 30.4 1.5 35.4 1.9 38.4 2.2 43.9 2.4 42.1 2.1 40.8 3.2 49.9 2.9 55.3 2.8
C23 5 2.6 18 3.7 15.4 4 14.5 4.5 21.9 2.8 16.1 4.7 26.3 4.5 31.9 3.7
# Remaining Off-Diagonal Moduli
C15 -2.3 0.3 -0.4 0.4 -0.2 0.4 0.1 0.5 1.2 0.4 2.3 1 3.2 0.6 5.1 0.6
C25 -7.8 0.7 -2.9 0.8 -5.1 1.1 -4.8 1.2 0.7 0.8 3.1 1.6 5.4 1 3.5 0.9
C35 7.5 0.6 4.6 0.8 7.2 1.1 6.9 1 2.5 0.8 2.2 1.5 1.7 0.9 0.5 0.9
C46 -7.2 0.1 -5.2 0.2 -4.8 0.2 -3.8 0.2 1.4 0.1 0.3 0.2 0.9 0.2 -2.2 0.1
C14 4.9 0.2 6.1 0.3 5.9 0.3 6.5 0.4 7.6 0.3 9.2 0.6 9 0.5 9.5 0.5
C16 -0.9 0.3 -0.6 0.4 -0.4 0.5 -0.8 0.5 -7.7 0.5 -10.1 1.4 -3 0.6 -10.8 0.7
C24 -4.4 0.6 -5.9 0.6 -7 0.6 -2.4 0.6 -2.9 0.5 0.9 1 2.1 0.9 7.5 0.6
C26 -6.4 0.9 -6.5 0.9 -6.8 1.2 -9.9 1.2 -6.8 1.1 -2.9 2.1 -9.9 1.3 -7.2 1.3
C34 -9.2 0.4 -2.9 0.5 2.2 0.7 -0.4 0.5 0.2 0.5 -0.9 1 1.7 0.9 6.6 0.6
C36 -9.4 0.6 -10.7 0.9 -9.8 0.9 -5.7 1 0.7 0.8 -0.3 1.2 -8.1 1.1 1.6 1
C45 -2.4 0.1 -1.3 0.1 -1.1 0.2 -1 0.2 0.2 0.1 0.7 0.2 0.8 0.1 3 0.1
C56 0.6 0.1 0.8 0.2 1.4 0.2 2.1 0.3 2.8 0.2 3.2 0.3 4.5 0.3 5.2 0.2
109 changes: 109 additions & 0 deletions contrib/anisotropic_plagioclase/plagioclase_data.py
View file
@@ -0,0 +1,109 @@
import pandas as pd
import numpy as np
from burnman.utils.unitcell import molar_volume_from_unit_cell_volume


def get_data():
d = pd.read_csv(
"data/plagioclase_stiffness_tensor_Brown_2016.dat",
comment="#",
header=0,
index_col=0,
sep="\\s+",
)

rows = list(d.index[1:])
d = np.array(d.to_numpy()[1:], dtype=float)
p_an = d[0, ::2] / 100.0 # fractions

CN = np.empty((len(p_an), 6, 6))
CN_err = np.empty((len(p_an), 6, 6))
for irow in range(1, 22):
name = rows[irow]
i, j = int(name[1]) - 1, int(name[2]) - 1
CN[:, i, j] = d[irow, ::2]
CN[:, j, i] = d[irow, ::2]
CN_err[:, i, j] = d[irow, 1::2]
CN_err[:, j, i] = d[irow, 1::2]

SN = np.linalg.inv(CN)
beta_RN = np.sum(SN[:, :3, :3], axis=(1, 2))
KRN = 1.0 / beta_RN
psiN = np.einsum("ijk, i->ijk", SN, KRN)
KRN_err = np.sqrt(np.sum(np.power(CN_err[:, :3, :3], 2.0), axis=(1, 2))) / 9.0
data = {
"CN": {
"P": p_an * 0.0 + 1.0e5,
"T": p_an * 0.0 + 298.15,
"p_an": p_an,
"psiN": psiN,
"SN": SN / 1.0e9,
"CN": CN * 1.0e9,
"CN_err": CN_err * 1.0e9,
"KRN": KRN * 1.0e9,
"KRN_err": KRN_err,
}
}

d = pd.read_csv(
"data/plagioclase_cell_tensor_Brown_2016.dat",
comment="#",
header=0,
index_col=0,
sep="\\s+",
).to_dict("list")
data["cell"] = {k: np.array(v) for k, v in d.items()}

# convert units
f_V = molar_volume_from_unit_cell_volume(1.0, 4.0)
f_ln = np.cbrt(f_V)
for prp in ["a", "a_err", "b", "b_err", "c", "c_err"]:
data["cell"][prp] = data["cell"][prp] * f_ln
for prp in ["V", "V_err"]:
data["cell"][prp] = data["cell"][prp] * f_V
for prp in ["rho", "rho_err"]:
data["cell"][prp] = data["cell"][prp] * 1.0e6

mask = data["cell"]["Z"] > 6.0
data["cell"]["c"][mask] = data["cell"]["c"][mask] / 2.0
data["cell"]["c_err"][mask] = data["cell"]["c_err"][mask] / 2.0
data["cell"]["V"][mask] = data["cell"]["V"][mask] / 2.0
data["cell"]["V_err"][mask] = data["cell"]["V_err"][mask] / 2.0

data["cell"]["P"] = data["cell"]["p_an"] * 0.0 + 1.0e5
data["cell"]["T"] = data["cell"]["p_an"] * 0.0 + 298.15

d = np.loadtxt("data/plagioclase_compressibilities_Brown_2016.dat")[:, 1:]

b = np.array([d[:, 2 * i + 1] for i in range(6)]) / 1.0e9
b_err = np.array([d[:, 2 * i + 2] for i in range(6)]) / 1.0e9
b = np.moveaxis(b, 1, 0)
b_err = np.moveaxis(b_err, 1, 0)

data["beta"] = {
"P": d[:, 0] * 0.0 + 1.0e5,
"T": d[:, 0] * 0.0 + 298.15,
"p_an": d[:, 0] / 100.0,
"b": b,
"b_err": b_err,
"bTR": np.sum(b[:, :3], axis=1),
"bTR_err": np.sqrt(np.sum(np.power(b_err[:, :3], 2.0), axis=1)),
}

return data


def get_C1_data():
data = get_data()
for dataset, d in data.items():
if dataset == "beta":
mask = d["p_an"] < 0.47
else:
mask = d["p_an"] < 0.50
data[dataset] = {k: v[mask] for k, v in d.items()}
return data


if __name__ == "__main__":
data = get_C1_data()
print(data)

0 comments on commit dd46dcb

Please sign in to comment.