From acb682e80d61c86d29875fa182d6a27a7f284e76 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Mon, 25 Mar 2024 23:54:00 +0000 Subject: [PATCH] anisotropic plagioclase --- contrib/anisotropic_plagioclase/Readme.md | 34 ++ .../plagioclase_cell_tensor_Brown_2016.dat | 14 + ...agioclase_compressibilities_Brown_2016.dat | 11 + ...lagioclase_stiffness_tensor_Brown_2016.dat | 26 ++ .../plagioclase_data.py | 109 ++++++ .../plagioclase_fit_eos.py | 213 +++++++++++ .../plagioclase_model.py | 192 ++++++++++ .../plagioclase_model_plots.py | 337 ++++++++++++++++++ .../plagioclase_parameters.py | 199 +++++++++++ 9 files changed, 1135 insertions(+) create mode 100644 contrib/anisotropic_plagioclase/Readme.md create mode 100644 contrib/anisotropic_plagioclase/data/plagioclase_cell_tensor_Brown_2016.dat create mode 100644 contrib/anisotropic_plagioclase/data/plagioclase_compressibilities_Brown_2016.dat create mode 100644 contrib/anisotropic_plagioclase/data/plagioclase_stiffness_tensor_Brown_2016.dat create mode 100644 contrib/anisotropic_plagioclase/plagioclase_data.py create mode 100644 contrib/anisotropic_plagioclase/plagioclase_fit_eos.py create mode 100644 contrib/anisotropic_plagioclase/plagioclase_model.py create mode 100644 contrib/anisotropic_plagioclase/plagioclase_model_plots.py create mode 100644 contrib/anisotropic_plagioclase/plagioclase_parameters.py diff --git a/contrib/anisotropic_plagioclase/Readme.md b/contrib/anisotropic_plagioclase/Readme.md new file mode 100644 index 00000000..2cb21e4c --- /dev/null +++ b/contrib/anisotropic_plagioclase/Readme.md @@ -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. diff --git a/contrib/anisotropic_plagioclase/data/plagioclase_cell_tensor_Brown_2016.dat b/contrib/anisotropic_plagioclase/data/plagioclase_cell_tensor_Brown_2016.dat new file mode 100644 index 00000000..1fe38a01 --- /dev/null +++ b/contrib/anisotropic_plagioclase/data/plagioclase_cell_tensor_Brown_2016.dat @@ -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 diff --git a/contrib/anisotropic_plagioclase/data/plagioclase_compressibilities_Brown_2016.dat b/contrib/anisotropic_plagioclase/data/plagioclase_compressibilities_Brown_2016.dat new file mode 100644 index 00000000..165a646e --- /dev/null +++ b/contrib/anisotropic_plagioclase/data/plagioclase_compressibilities_Brown_2016.dat @@ -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 diff --git a/contrib/anisotropic_plagioclase/data/plagioclase_stiffness_tensor_Brown_2016.dat b/contrib/anisotropic_plagioclase/data/plagioclase_stiffness_tensor_Brown_2016.dat new file mode 100644 index 00000000..69ad7232 --- /dev/null +++ b/contrib/anisotropic_plagioclase/data/plagioclase_stiffness_tensor_Brown_2016.dat @@ -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 diff --git a/contrib/anisotropic_plagioclase/plagioclase_data.py b/contrib/anisotropic_plagioclase/plagioclase_data.py new file mode 100644 index 00000000..4b2c84b3 --- /dev/null +++ b/contrib/anisotropic_plagioclase/plagioclase_data.py @@ -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) diff --git a/contrib/anisotropic_plagioclase/plagioclase_fit_eos.py b/contrib/anisotropic_plagioclase/plagioclase_fit_eos.py new file mode 100644 index 00000000..6cf74e56 --- /dev/null +++ b/contrib/anisotropic_plagioclase/plagioclase_fit_eos.py @@ -0,0 +1,213 @@ +import numpy as np +from plagioclase_model import make_scalar_model, make_anisotropic_model +from plagioclase_data import get_C1_data +from scipy.optimize import minimize +from plagioclase_parameters import scalar_args, cell_args, elastic_args +from scipy.optimize import differential_evolution + +data = get_C1_data() + +min_chisqr = [1.0e18] +iteration = [0] + + +def misfit_scalar(args): + # For scalar EoS, fit isothermal compressibilities and volumes + ss = make_scalar_model(args) + + molar_fractions = np.array([data["cell"]["p_an"], 1.0 - data["cell"]["p_an"]]).T + pressures = data["cell"]["P"] + temperatures = data["cell"]["T"] + + V_obs = ss.evaluate(["V"], pressures, temperatures, molar_fractions)[0] + + chisqr = np.sum( + np.power((V_obs - data["cell"]["V"]) / (data["cell"]["V_err"]), 2.0) + ) + + molar_fractions = np.array([data["beta"]["p_an"], 1.0 - data["beta"]["p_an"]]).T + pressures = data["beta"]["P"] + temperatures = data["beta"]["T"] + + beta_model = ss.evaluate( + ["isothermal_compressibility_reuss"], pressures, temperatures, molar_fractions + )[0] + chisqr += np.sum( + np.power((beta_model - data["beta"]["bTR"]) / (data["beta"]["bTR_err"]), 2.0) + ) + + # add a chisqr for extreme nonlinear mixing behaviour in K_T + # chisqr += np.power(args[5]/4., 2.) + + if chisqr < min_chisqr[0]: + min_chisqr[0] = chisqr + print(repr(args)) + print(chisqr) + print() + return chisqr + + +def misfit_cell(args, scalar_args, elastic_args): + # For cell EoS, fit cell parameters and isothermal compressibilities + + cell_args = args + ss = make_anisotropic_model(scalar_args, cell_args, elastic_args) + + # cell parameters + pressures = data["cell"]["P"] + temperatures = data["cell"]["T"] + molar_fractions = np.array([data["cell"]["p_an"], 1.0 - data["cell"]["p_an"]]).T + cell_model = ss.evaluate( + ["cell_parameters"], pressures, temperatures, molar_fractions + )[0] + + chisqr = 0.0 + f_cell = 0.1 # cell and compressibility uncertainties are poorly scaled + for i, prm in enumerate(["a", "b", "c", "alpha", "beta", "gamma"]): + chi = (cell_model[:, i] - data["cell"][prm]) / data["cell"][prm + "_err"] + chisqr += f_cell * np.sum(np.power(chi, 2.0)) + + # isothermal compressibilities + if False: + pressures = data["beta"]["P"] + temperatures = data["beta"]["T"] + molar_fractions = np.array([data["beta"]["p_an"], 1.0 - data["beta"]["p_an"]]).T + beta_model = ss.evaluate( + ["isothermal_compressibility_tensor"], + pressures, + temperatures, + molar_fractions, + )[0] + + f_beta = 10.0 # cell and compressibility uncertainties are poorly scaled + for idx, (i, j) in enumerate([(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)]): + beta_model_i = beta_model[:, i, j] + chi = (beta_model_i - data["beta"]["b"][:, idx]) / data["beta"]["b_err"][ + :, idx + ] + chisqr += f_beta * np.sum(np.power(chi, 2.0)) + + if False: + iteration[0] = iteration[0] + 1 + if chisqr < min_chisqr[0]: + min_chisqr[0] = chisqr + print(repr(args)) + print(chisqr) + else: + print(f"{iteration[0]}: {chisqr} \r") + + return chisqr + + +def misfit_elastic(args, scalar_args, cell_args): + elastic_args = args + ss = make_anisotropic_model(scalar_args, cell_args, elastic_args) + + pressures = data["CN"]["P"] + temperatures = data["CN"]["T"] + molar_fractions = np.array([data["CN"]["p_an"], 1.0 - data["CN"]["p_an"]]).T + CN_model = ss.evaluate( + ["isentropic_stiffness_tensor"], pressures, temperatures, molar_fractions + )[0] + + chi = (CN_model - data["CN"]["CN"]) / data["CN"]["CN_err"] + chisqr = np.sum(np.tri(6, 6) * np.power(chi, 2.0)) + + if False: + iteration[0] = iteration[0] + 1 + if chisqr < min_chisqr[0]: + min_chisqr[0] = chisqr + print(repr(args)) + print(chisqr) + else: + print(f"\r{iteration[0]}: {chisqr}", end="", flush=True) + + return chisqr + + +def misfit_cell_and_elastic(args, scalar_args): + cell_args = args[:20] + elastic_args = args[20:] + chisqr = misfit_cell(cell_args, scalar_args, elastic_args) + chisqr += misfit_elastic(elastic_args, scalar_args, cell_args) + + iteration[0] = iteration[0] + 1 + if chisqr < min_chisqr[0]: + min_chisqr[0] = chisqr + print(repr(args[:20])) + print(repr(args[20:])) + print(chisqr) + else: + print(f"\r{iteration[0]}: {chisqr}", end="", flush=True) + + return chisqr + + +def misfit_elastic_2(args, scalar_args, cell_args): + elastic_args = args + ss = make_anisotropic_model(scalar_args, cell_args, elastic_args) + + pressures = data["CN"]["P"] + temperatures = data["CN"]["T"] + molar_fractions = np.array([data["CN"]["p_an"], 1.0 - data["CN"]["p_an"]]).T + SN_model, KNR_model = ss.evaluate( + ["isentropic_compliance_tensor", "isentropic_bulk_modulus_reuss"], + pressures, + temperatures, + molar_fractions, + ) + phi_model = np.einsum("ijk, i->ijk", SN_model, KNR_model) + + chi = phi_model - data["CN"]["phiN"] + chisqr = np.sum(np.power(chi, 2.0)) + + iteration[0] = iteration[0] + 1 + if chisqr < min_chisqr[0]: + min_chisqr[0] = chisqr + print(repr(args)) + print(chisqr) + else: + print(f"\r{iteration[0]}: {chisqr}", end="", flush=True) + + return chisqr + + +if False: + args = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) + sol = minimize(misfit_scalar, args) + exit() + +if False: + sol = minimize( + misfit_cell, + cell_args, + args=(scalar_args, elastic_args), + ) + +if False: + bounds = [(-20, 20) for i in range(30)] + sol = minimize( + misfit_elastic, + elastic_args, + args=(scalar_args, cell_args), + ) # bounds=bounds) + +if False: + bounds = [(-10, 10) for i in range(30)] + sol = differential_evolution( + misfit_elastic, bounds=bounds, x0=elastic_args, args=(scalar_args, cell_args) + ) + +if True: + ce_args = np.concatenate((cell_args, elastic_args)) + sol = minimize(misfit_cell_and_elastic, ce_args, args=(scalar_args)) + + sol = minimize( + misfit_cell_and_elastic, sol.x, args=(scalar_args), method="Nelder-Mead" + ) + + sol = minimize( + misfit_cell_and_elastic, sol.x, args=(scalar_args), method="Nelder-Mead" + ) + + sol = minimize(misfit_cell_and_elastic, sol.x, args=(scalar_args)) diff --git a/contrib/anisotropic_plagioclase/plagioclase_model.py b/contrib/anisotropic_plagioclase/plagioclase_model.py new file mode 100644 index 00000000..82b5385d --- /dev/null +++ b/contrib/anisotropic_plagioclase/plagioclase_model.py @@ -0,0 +1,192 @@ +from burnman import Solution +from burnman.minerals import SLB_2022 +from burnman.classes.solutionmodel import PolynomialSolution +from burnman import AnisotropicMineral, AnisotropicSolution +import numpy as np +from copy import deepcopy +from burnman.utils.unitcell import cell_parameters_to_vectors + + +def make_scalar_model(args): + dV_ab, dV_an, dV_aban = args[:3] / 1.0e6 + dK_ab, dK_an, dK_aban = args[3:6] * 1.0e9 + + ab_scalar = SLB_2022.albite() + an_scalar = SLB_2022.anorthite() + + aban_linear = SLB_2022.anorthite() + aban_linear.params["V_0"] = 0.5 * ( + ab_scalar.params["V_0"] + an_scalar.params["V_0"] + ) + aban_linear.params["K_0"] = 0.5 * ( + ab_scalar.params["K_0"] + an_scalar.params["K_0"] + ) + + aban_scalar = deepcopy(aban_linear) + + ab_scalar.params["V_0"] = ab_scalar.params["V_0"] + dV_ab + an_scalar.params["V_0"] = an_scalar.params["V_0"] + dV_an + aban_scalar.params["V_0"] = aban_scalar.params["V_0"] + dV_aban + + ab_scalar.params["K_0"] = ab_scalar.params["K_0"] + dK_ab + an_scalar.params["K_0"] = an_scalar.params["K_0"] + dK_an + aban_scalar.params["K_0"] = aban_scalar.params["K_0"] + dK_aban + + class plagioclase_scalar(Solution): + def __init__(self, molar_fractions=None): + self.name = "plagioclase (NCAS)" + self.solution_model = PolynomialSolution( + endmembers=[ + [an_scalar, "[Ca][Al]2Si2O8"], + [ab_scalar, "[Na][Al1/2Si1/2]2Si2O8"], + ], + interaction_endmembers=[aban_scalar, aban_linear], + endmember_coefficients_and_interactions=[[4.0, -4.0, 0, 1, 1, 1]], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + return plagioclase_scalar() + + +def make_anisotropic_model(scalar_args, cell_args, elastic_args): + + ss = make_scalar_model(scalar_args) + an_scalar = ss.endmembers[0][0] + ab_scalar = ss.endmembers[1][0] + aban_scalar = ss.solution_model.interaction_endmembers[0] + aban_linear = ss.solution_model.interaction_endmembers[1] + + b_over_a_ab, c_over_a_ab = cell_args[:2] + alpha_ab, beta_ab, gamma_ab = cell_args[2:5] + ab_cell_parameters = np.array( + [1.0, b_over_a_ab, c_over_a_ab, alpha_ab, beta_ab, gamma_ab] + ) + + frame_convention = [1, 2, 0] + + M = cell_parameters_to_vectors(ab_cell_parameters, frame_convention) + f = np.cbrt(ab_scalar.params["V_0"] / np.linalg.det(M)) + ab_cell_parameters[:3] = ab_cell_parameters[:3] * f + + b_over_a_an, c_over_a_an = cell_args[5:7] + alpha_an, beta_an, gamma_an = cell_args[7:10] + an_cell_parameters = np.array( + [1.0, b_over_a_an, c_over_a_an, alpha_an, beta_an, gamma_an] + ) + M = cell_parameters_to_vectors(an_cell_parameters, frame_convention) + f = np.cbrt(an_scalar.params["V_0"] / np.linalg.det(M)) + an_cell_parameters[:3] = an_cell_parameters[:3] * f + + an_a = np.zeros((6, 6)) + ab_a = np.zeros((6, 6)) + + # the following are all the elastic arguments + # 15x2 = 30 in total + # corresponding to all the block off-diagonals + idx = 0 + # 3*2 = 6; 44, 55, 66 + for i in range(3, 6): + an_a[i, i] = elastic_args[idx] + ab_a[i, i] = elastic_args[idx + 1] + idx = idx + 2 + # 6*2 = 12; 12, 13, 23, 15, 16, 26 + for j_shift in [0, 3]: + for i in range(2): + for j in range(i + 1, 3): + an_a[i, j + j_shift] = elastic_args[idx] + ab_a[i, j + j_shift] = elastic_args[idx + 1] + an_a[j + j_shift, i] = elastic_args[idx] + ab_a[j + j_shift, i] = elastic_args[idx + 1] + idx = idx + 2 + + # 6*2 = 12; 24, 34, 35, 45, 46, 56 + for i_shift in [0, 3]: + for i in range(1, 3): + for j in range(3, i + 3): + an_a[i + i_shift, j] = elastic_args[idx] + ab_a[i + i_shift, j] = elastic_args[idx + 1] + an_a[j, i + i_shift] = elastic_args[idx] + ab_a[j, i + i_shift] = elastic_args[idx + 1] + idx = idx + 2 + + # the following is another 5 x 2 = 10 cell arguments + # corresponding to all the block diagonals + # bringing the total to idx + 10 = 20 + an_b = np.sum(an_a[:3, :3], axis=0) + an_c = np.sum(an_a[:3, 3:], axis=0) + ab_b = np.sum(ab_a[:3, :3], axis=0) + ab_c = np.sum(ab_a[:3, 3:], axis=0) + + idx = 10 + # 2x2 = 4; 22, 33 + for i in range(1, 3): + j = i + an_a[i, i] = cell_args[idx] - an_b[i] + ab_a[i, i] = cell_args[idx + 1] - ab_b[i] + idx = idx + 2 + # 3x2 = 6; 14, 25, 36 + for i in range(3): + j = i + 3 + an_a[i, j] = cell_args[idx] - an_c[i] + ab_a[i, j] = cell_args[idx + 1] - ab_c[i] + an_a[j, i] = cell_args[idx] - an_c[i] + ab_a[j, i] = cell_args[idx + 1] - ab_c[i] + idx = idx + 2 + # Finally, 11 + an_a[0, 0] = 1.0 - np.sum(an_a[:3, :3]) + ab_a[0, 0] = 1.0 - np.sum(ab_a[:3, :3]) + + # print(an_a) + # print(ab_a) + # exit() + + an_anisotropic_parameters = {"a": an_a} + ab_anisotropic_parameters = {"a": ab_a} + + def psi_func(f, Pth, params): + dPsidf = params["a"] + Psi = params["a"] * f + dPsidPth = np.zeros((6, 6)) + return (Psi, dPsidf, dPsidPth) + + an_aniso = AnisotropicMineral( + an_scalar, + an_cell_parameters, + an_anisotropic_parameters, + psi_function=psi_func, + frame_convention=frame_convention, + orthotropic=False, + ) + + ab_aniso = AnisotropicMineral( + ab_scalar, + ab_cell_parameters, + ab_anisotropic_parameters, + psi_function=psi_func, + frame_convention=frame_convention, + orthotropic=False, + ) + + # An ideal mixing model with two endmembers + def ideal_fn(lnV, Pth, X, params): + z = np.zeros((6, 6)) + f = np.zeros((3, 3, 2)) + return (z, z, z, f) + + prm = {} + + plagioclase_anisotropic = AnisotropicSolution( + name="plagioclase (C1)", + solution_model=PolynomialSolution( + endmembers=[ + [an_aniso, "[Ca][Al]2Si2O8"], + [ab_aniso, "[Na][Al1/2Si1/2]2Si2O8"], + ], + interaction_endmembers=[aban_scalar, aban_linear], + endmember_coefficients_and_interactions=[[4.0, -4.0, 0, 1, 1, 1]], + ), + psi_excess_function=ideal_fn, + anisotropic_parameters=prm, + ) + + return plagioclase_anisotropic diff --git a/contrib/anisotropic_plagioclase/plagioclase_model_plots.py b/contrib/anisotropic_plagioclase/plagioclase_model_plots.py new file mode 100644 index 00000000..e6b372b6 --- /dev/null +++ b/contrib/anisotropic_plagioclase/plagioclase_model_plots.py @@ -0,0 +1,337 @@ +import numpy as np +import matplotlib.pyplot as plt +from plagioclase_parameters import scalar_args, cell_args, elastic_args +from plagioclase_model import make_anisotropic_model +from plagioclase_data import get_data + +ss = make_anisotropic_model(scalar_args, cell_args, elastic_args) + +# Model data +p_ans = np.linspace(0.0, 1.0, 101) +molar_fractions = np.array([p_ans, 1.0 - p_ans]).T +pressures = p_ans * 0.0 + 1.0e5 +temperatures = p_ans * 0.0 + 300.0 +prps = ss.evaluate( + [ + "molar_volume", + "isothermal_bulk_modulus_reuss", + "isothermal_stiffness_tensor", + "isentropic_stiffness_tensor", + "isothermal_compliance_tensor", + "isothermal_compressibility_tensor", + "cell_parameters", + ], + pressures, + temperatures, + molar_fractions, +) + +labels = ["V", "KTR", "CT", "CN", "ST", "betaT", "cell"] +prps = {labels[i]: prps[i] for i in range(7)} +prps["psi"] = np.einsum("ijk, i->ijk", prps["ST"], prps["KTR"]) + +d = get_data() + + +# Scalar properties +fig = plt.figure(figsize=(8, 4)) +ax = [fig.add_subplot(1, 2, i) for i in range(1, 3)] + + +mask2 = p_ans < 0.5 + +ln = ax[0].plot(p_ans[mask2], prps["V"][mask2] * 1.0e6) +ax[0].plot( + p_ans[~mask2], prps["V"][~mask2] * 1.0e6, color=ln[0].get_color(), linestyle=":" +) + +mask = d["cell"]["p_an"] < 0.5 +ax[0].scatter( + d["cell"]["p_an"][mask], d["cell"]["V"][mask] * 1.0e6, color=ln[0].get_color() +) +ax[0].errorbar( + d["cell"]["p_an"][mask], + d["cell"]["V"][mask] * 1.0e6, + d["cell"]["V_err"][mask] * 1.0e6, + color=ln[0].get_color(), + ls="none", +) +ax[0].scatter( + d["cell"]["p_an"][~mask], + d["cell"]["V"][~mask] * 1.0e6, + color=ln[0].get_color(), + alpha=0.5, +) +ax[0].errorbar( + d["cell"]["p_an"][~mask], + d["cell"]["V"][~mask] * 1.0e6, + d["cell"]["V_err"][~mask] * 1.0e6, + color=ln[0].get_color(), + alpha=0.5, + ls="none", +) + + +ax[1].plot(p_ans[mask2], prps["KTR"][mask2] / 1.0e9) +ax[1].plot( + p_ans[~mask2], prps["KTR"][~mask2] / 1.0e9, color=ln[0].get_color(), linestyle=":" +) + +mask = d["beta"]["p_an"] < 0.47 +ax[1].scatter( + d["beta"]["p_an"][mask], + 1.0 / d["beta"]["bTR"][mask] / 1.0e9, + color=ln[0].get_color(), +) +ax[1].errorbar( + d["beta"]["p_an"][mask], + 1.0 / d["beta"]["bTR"][mask] / 1.0e9, + (d["beta"]["bTR_err"] / d["beta"]["bTR"] / d["beta"]["bTR"])[mask] / 1.0e9, + color=ln[0].get_color(), + ls="none", +) +ax[1].scatter( + d["beta"]["p_an"][~mask], + 1.0 / d["beta"]["bTR"][~mask] / 1.0e9, + color=ln[0].get_color(), + alpha=0.5, +) +ax[1].errorbar( + d["beta"]["p_an"][~mask], + 1.0 / d["beta"]["bTR"][~mask] / 1.0e9, + (d["beta"]["bTR_err"] / d["beta"]["bTR"] / d["beta"]["bTR"])[~mask] / 1.0e9, + color=ln[0].get_color(), + alpha=0.5, + ls="none", +) + +for i in range(2): + ax[i].set_xlabel("$p_{an}$") + +ax[0].set_ylabel("$V$ (m$^3$/mol)") +ax[1].set_ylabel("$K_{\\text{TR}}$ (GPa)") + +fig.set_tight_layout(True) +fig.savefig("plag_V_KT.pdf") + +# Elastic stiffness data +ijs = np.array( + [ + [1, 1], + [2, 2], + [3, 3], + [4, 4], + [5, 5], + [6, 6], + [1, 2], + [1, 3], + [2, 3], + [1, 5], + [2, 5], + [3, 5], + [4, 6], + [1, 4], + [1, 6], + [2, 4], + [2, 6], + [3, 4], + [3, 6], + [4, 5], + [5, 6], + ] +) + +inds = ijs - 1 +mask = d["CN"]["p_an"] < 0.5 +mask2 = p_ans < 0.5 + + +fig = plt.figure(figsize=(10, 8)) +ax = [fig.add_subplot(3, 3, i) for i in range(1, 8)] +for irow, (i, j) in enumerate(inds): + name = f"$C_{{{i+1}{j+1}}}$" + axi = int((irow - (irow) % 3) / 3) + ln = ax[axi].plot(p_ans[mask2], prps["CN"][mask2, i, j] / 1.0e9) + ax[axi].plot( + p_ans[~mask2], + prps["CN"][~mask2, i, j] / 1.0e9, + linestyle=":", + color=ln[0].get_color(), + ) + ax[axi].scatter( + d["CN"]["p_an"][mask], + d["CN"]["CN"][mask, i, j] / 1.0e9, + label=name, + color=ln[0].get_color(), + ) + ax[axi].errorbar( + d["CN"]["p_an"][mask], + d["CN"]["CN"][mask, i, j] / 1.0e9, + d["CN"]["CN_err"][mask, i, j] / 1.0e9, + ls="none", + color=ln[0].get_color(), + ) + ax[axi].scatter( + d["CN"]["p_an"][~mask], + d["CN"]["CN"][~mask, i, j] / 1.0e9, + color=ln[0].get_color(), + alpha=0.5, + ) + ax[axi].errorbar( + d["CN"]["p_an"][~mask], + d["CN"]["CN"][~mask, i, j] / 1.0e9, + d["CN"]["CN_err"][~mask, i, j] / 1.0e9, + ls="none", + color=ln[0].get_color(), + alpha=0.5, + ) + +for i in range(7): + ax[i].legend() + ax[i].set_xlabel("$p_{an}$") + ax[i].set_ylabel("modulus (GPa)") + +fig.set_tight_layout(True) +fig.savefig("plag_stiffnesses.pdf") + +# psi +fig = plt.figure(figsize=(10, 8)) +ax = [fig.add_subplot(3, 3, i) for i in range(1, 8)] +for irow, (i, j) in enumerate(inds): + name = f"$d\\Psi_{{{i+1}{j+1}}} / df$" + axi = int((irow - (irow) % 3) / 3) + ln = ax[axi].plot(p_ans[mask2], prps["psi"][mask2, i, j]) + ax[axi].plot( + p_ans[~mask2], prps["psi"][~mask2, i, j], linestyle=":", color=ln[0].get_color() + ) + ax[axi].scatter( + d["CN"]["p_an"][mask], + d["CN"]["psiN"][mask, i, j], + label=name, + color=ln[0].get_color(), + ) + ax[axi].scatter( + d["CN"]["p_an"][~mask], + d["CN"]["psiN"][~mask, i, j], + color=ln[0].get_color(), + alpha=0.5, + ) + +for i in range(7): + ax[i].legend() + ax[i].set_xlabel("$p_{an}$") + ax[i].set_ylabel("$d\\Psi / df$") + +fig.set_tight_layout(True) +fig.savefig("plag_psi.pdf") + +# Cell parameters +fig = plt.figure(figsize=(10, 6)) +ax = [fig.add_subplot(2, 3, i) for i in range(1, 7)] + +labels = ["a", "b", "c", "alpha", "beta", "gamma"] + +inds = ijs - 1 +mask = d["cell"]["p_an"] < 0.5 +mask2 = p_ans < 0.5 +for i, label in enumerate(labels): + + ln = ax[i].plot(p_ans[mask2], prps["cell"][mask2, i]) + ax[i].plot( + p_ans[~mask2], prps["cell"][~mask2, i], linestyle=":", color=ln[0].get_color() + ) + ax[i].scatter( + d["cell"]["p_an"][mask], d["cell"][label][mask], color=ln[0].get_color() + ) + ax[i].errorbar( + d["cell"]["p_an"][mask], + d["cell"][label][mask], + d["cell"][label + "_err"][mask], + ls="none", + color=ln[0].get_color(), + ) + ax[i].scatter( + d["cell"]["p_an"][~mask], + d["cell"][label][~mask], + color=ln[0].get_color(), + alpha=0.5, + ) + ax[i].errorbar( + d["cell"]["p_an"][~mask], + d["cell"][label][~mask], + d["cell"][label + "_err"][~mask], + ls="none", + color=ln[0].get_color(), + alpha=0.5, + ) + + if i < 3: + ax[i].set_ylabel(f"${label}$ (m)") + else: + ax[i].set_ylabel(f"$\\{label}$ ($^{{\\circ}}$)") + ax[i].set_xlabel("$p_{an}$") + +fig.set_tight_layout(True) +fig.savefig("plag_cell_parameters.pdf") + +# Compressibilities +fig = plt.figure(figsize=(10, 6)) +ax = [fig.add_subplot(2, 3, i) for i in range(1, 7)] + +ijs = np.array( + [ + [1, 1], + [2, 2], + [3, 3], + [2, 3], + [1, 3], + [1, 2], + ] +) + +inds = ijs - 1 + +mask = d["beta"]["p_an"] < 0.5 +mask2 = p_ans < 0.5 +for axi, (i, j) in enumerate(inds): + + ln = ax[axi].plot(p_ans[mask2], prps["betaT"][mask2, i, j] * 1.0e9) + ax[axi].plot( + p_ans[~mask2], + prps["betaT"][~mask2, i, j] * 1.0e9, + linestyle=":", + color=ln[0].get_color(), + ) + ax[axi].scatter( + d["beta"]["p_an"][mask], + d["beta"]["b"][mask, axi] * 1.0e9, + color=ln[0].get_color(), + ) + ax[axi].errorbar( + d["beta"]["p_an"][mask], + d["beta"]["b"][mask, axi] * 1.0e9, + d["beta"]["b_err"][mask, axi] * 1.0e9, + ls="none", + color=ln[0].get_color(), + ) + ax[axi].scatter( + d["beta"]["p_an"][~mask], + d["beta"]["b"][~mask, axi] * 1.0e9, + color=ln[0].get_color(), + alpha=0.5, + ) + ax[axi].errorbar( + d["beta"]["p_an"][~mask], + d["beta"]["b"][~mask, axi] * 1.0e9, + d["beta"]["b_err"][~mask, axi] * 1.0e9, + ls="none", + color=ln[0].get_color(), + alpha=0.5, + ) + + ax[axi].set_ylabel(f"$\\beta_{{{axi}}}$ (m)") + ax[axi].set_xlabel("$p_{an}$") + +fig.set_tight_layout(True) +fig.savefig("plag_compressibilities.pdf") +plt.show() diff --git a/contrib/anisotropic_plagioclase/plagioclase_parameters.py b/contrib/anisotropic_plagioclase/plagioclase_parameters.py new file mode 100644 index 00000000..c23f736d --- /dev/null +++ b/contrib/anisotropic_plagioclase/plagioclase_parameters.py @@ -0,0 +1,199 @@ +import numpy as np + +scalar_args = np.array( + [-0.48217651, 0.5647768, 0.18958683, -4.53426789, 4.35928, 7.19206864] +) + +cell_args = np.array( + [ + 1.57145311e00, + 8.79492227e-01, + 9.42490900e01, + 1.16600835e02, + 8.77870251e01, + 1.57522940e00, + 8.59890806e-01, + 9.30305623e01, + 1.16096344e02, + 9.22130388e01, + 7.45685978e-02, + 2.55743752e-01, + 5.95265343e-01, + 1.38172513e-01, + -1.06464568e00, + 2.32867395e-01, + -3.85637532e-01, + 1.14771229e-01, + 1.19734522e00, + -6.36213240e-02, + ] +) + +cell_args = np.array( + [ + 1.57145311e00, + 8.79492227e-01, + 9.42490900e01, + 1.16600835e02, + 8.77870251e01, + 1.57522940e00, + 8.59890806e-01, + 9.30305623e01, + 1.16096344e02, + 9.22130388e01, + 0.27, + 0.1851, + 0.421, + 0.193, + -0.2 * 2, + 0.02 * 2.0, + -0.125 * 2.0, + 0.05 * 2.0, + 0.15 * 2.0, + 0.11 * 2.0, + ] +) + +cell_args = [ + 1.57138375e00, + 8.80353680e-01, + 9.42621730e01, + 1.16625876e02, + 8.77842388e01, + 1.57499235e00, + 8.58223150e-01, + 9.29880407e01, + 1.16105350e02, + 9.22823088e01, + 0.27, + 0.185, + 0.421, + 0.193, + -6.06397972e-01, + 1.17616351e-01, + -2.97717520e-01, + 3.52587083e-02, + 4.91428406e-01, + 1.62254178e-01, +] + +cell_args = [ + 1.57135650e00, + 8.79690034e-01, + 9.42681953e01, + 1.16653731e02, + 8.77456286e01, + 1.57727574e00, + 8.58786646e-01, + 9.29321860e01, + 1.16101263e02, + 9.24015066e01, + 0.27, + 0.185, + 0.421, + 0.193, + -6.69965921e-01, + 1.54062702e-01, + -3.75783217e-01, + 6.48340839e-02, + 5.84238460e-01, + 5.22397885e-03, +] + +# 0.122 +elastic_args = np.array( + [ + 4.33909620e00, + 2.54092660e00, + 2.42175356e00, + 2.07815062e00, + 2.56071749e00, + 1.78316916e00, + -2.99841390e-01, + -1.86664570e-01, + -2.41161330e-01, + -1.80809767e-01, + 3.68197415e-03, + 2.85613898e-02, + -1.17713816e-01, + 6.87778295e-02, + -7.39152187e-02, + -9.50888511e-02, + 2.00337451e-01, + 1.31234752e-01, + 1.36714207e-01, + 1.63787559e-01, + -1.91896776e-01, + 2.09237554e-01, + -1.63083822e-01, + -6.72865617e-02, + -3.94445929e-02, + 1.74636863e-01, + 1.95199461e-01, + 6.16363958e-01, + -2.45633713e-01, + -7.66198811e-04, + ] +) + +# 1624 combined +cell_args = np.array( + [ + 1.57151918e00, + 8.79431200e-01, + 9.42442794e01, + 1.16593320e02, + 8.77955924e01, + 1.57782743e00, + 8.59025416e-01, + 9.26631785e01, + 1.16006800e02, + 9.27426526e01, + 2.06560097e-01, + 2.04196306e-01, + 4.54456435e-01, + 1.84459426e-01, + -2.36438516e-01, + 3.12910093e-02, + -2.09187316e-01, + 6.30106933e-02, + 1.93416913e-01, + 9.07804895e-02, + ] +) + + +elastic_args = np.array( + [ + 4.27237738, + 2.59013263, + 2.41180097, + 2.1290894, + 2.56196346, + 1.82347664, + -0.31242721, + -0.19103324, + -0.18527108, + -0.20792989, + -0.02411557, + 0.03255825, + -0.06669763, + 0.05813644, + -0.09332946, + -0.14800036, + 0.16297612, + 0.09251345, + 0.19836256, + 0.16537785, + -0.12707124, + 0.21482192, + -0.11676085, + -0.08675306, + -0.06509725, + 0.1760466, + 0.19217078, + 0.62808897, + -0.24125105, + -0.007342, + ] +)