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.

Sign up for GitHub

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

added tabulation for plag contribution #590

Merged
merged 1 commit into from
May 19, 2024
Merged

added tabulation for plag contribution #590

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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: 4 additions & 0 deletions contrib/anisotropic_plagioclase/plagioclase_data.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,10 @@ def get_data():
b = np.moveaxis(b, 1, 0)
b_err = np.moveaxis(b_err, 1, 0)

# Brown reports the sum in Voigt form (Section 3.1)
b[:, 3:] = b[:, 3:] / 2
b_err[:, 3:] = b_err[:, 3:] / 2

data["beta"] = {
"P": d[:, 0] * 0.0 + 1.0e5,
"T": d[:, 0] * 0.0 + 298.15,
Expand Down
49 changes: 44 additions & 5 deletions contrib/anisotropic_plagioclase/plagioclase_model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
import numpy as np
from copy import deepcopy
from burnman.utils.unitcell import cell_parameters_to_vectors
from tabulate import tabulate


def make_scalar_model(args):
Expand All @@ -14,24 +15,46 @@ def make_scalar_model(args):
ab_scalar = SLB_2022.albite()
an_scalar = SLB_2022.anorthite()

ab_scalar.params["V_0"] = ab_scalar.params["V_0"] + dV_ab
an_scalar.params["V_0"] = an_scalar.params["V_0"] + dV_an
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_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_linear.params["Kprime_0"] = 0.5 * (
ab_scalar.params["Kprime_0"] + an_scalar.params["Kprime_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

table = [
["", "ab", "an", "an$_{50}$ (1)", "an$_{50}$ (2)"],
[
"V_0 (cm$^3$/mol)",
ab_scalar.params["V_0"] * 1.0e6,
an_scalar.params["V_0"] * 1.0e6,
aban_linear.params["V_0"] * 1.0e6,
aban_scalar.params["V_0"] * 1.0e6,
],
[
"K_0 (GPa)",
ab_scalar.params["K_0"] / 1.0e9,
an_scalar.params["K_0"] / 1.0e9,
aban_linear.params["K_0"] / 1.0e9,
aban_scalar.params["K_0"] / 1.0e9,
],
]
print(tabulate(table, headers="firstrow", tablefmt="latex_raw", floatfmt=".6e"))

class plagioclase_scalar(Solution):
def __init__(self, molar_fractions=None):
self.name = "plagioclase (NCAS)"
Expand Down Expand Up @@ -77,6 +100,16 @@ def make_anisotropic_model(scalar_args, cell_args, elastic_args):
f = np.cbrt(an_scalar.params["V_0"] / np.linalg.det(M))
an_cell_parameters[:3] = an_cell_parameters[:3] * f

table = [
["", "$a$", "$b$", "$c$", "$\\alpha$", "$\\beta$", "$\\gamma$"],
["ab"],
["an"],
]
table[1].extend(ab_cell_parameters)
table[2].extend(an_cell_parameters)

print(tabulate(table, headers="firstrow", tablefmt="latex_raw", floatfmt=".6e"))

an_a = np.zeros((6, 6))
ab_a = np.zeros((6, 6))

Expand Down Expand Up @@ -136,6 +169,12 @@ def make_anisotropic_model(scalar_args, cell_args, elastic_args):
an_a[0, 0] = 1.0 - np.sum(an_a[:3, :3])
ab_a[0, 0] = 1.0 - np.sum(ab_a[:3, :3])

table = ab_a
print(tabulate(table, tablefmt="latex_raw", floatfmt=".5f"))

table = an_a
print(tabulate(table, tablefmt="latex_raw", floatfmt=".5f"))

# print(an_a)
# print(ab_a)
# exit()
Expand Down
38 changes: 31 additions & 7 deletions contrib/anisotropic_plagioclase/plagioclase_model_plots.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,9 @@
from plagioclase_parameters import scalar_args, cell_args, elastic_args
from plagioclase_model import make_anisotropic_model
from plagioclase_data import get_data
from burnman.minerals.SLB_2022 import plagioclase

ss_SLB = plagioclase()

ss = make_anisotropic_model(scalar_args, cell_args, elastic_args)

Expand All @@ -26,10 +29,20 @@
molar_fractions,
)

prps_SLB = ss_SLB.evaluate(
["molar_volume", "isothermal_bulk_modulus_reuss"],
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"])

prps_SLB = {labels[i]: prps_SLB[i] for i in range(2)}

d = get_data()


Expand All @@ -40,7 +53,15 @@

mask2 = p_ans < 0.5

ln = ax[0].plot(p_ans[mask2], prps["V"][mask2] * 1.0e6)
ln = ax[0].plot(p_ans[mask2], prps["V"][mask2] * 1.0e6, label="C$\\bar{1}$, this study")
ln = ax[0].plot(
p_ans,
prps_SLB["V"] * 1.0e6,
linestyle="--",
color=ln[0].get_color(),
label="SLB2022",
)

ax[0].plot(
p_ans[~mask2], prps["V"][~mask2] * 1.0e6, color=ln[0].get_color(), linestyle=":"
)
Expand Down Expand Up @@ -72,7 +93,8 @@
)


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())
ax[1].plot(p_ans, prps_SLB["KTR"] / 1.0e9, linestyle="--", color=ln[0].get_color())
ax[1].plot(
p_ans[~mask2], prps["KTR"][~mask2] / 1.0e9, color=ln[0].get_color(), linestyle=":"
)
Expand Down Expand Up @@ -108,9 +130,11 @@
for i in range(2):
ax[i].set_xlabel("$p_{an}$")

ax[0].set_ylabel("$V$ (m$^3$/mol)")
ax[0].set_ylabel("$V$ (cm$^3$/mol)")
ax[1].set_ylabel("$K_{\\text{TR}}$ (GPa)")

ax[0].legend()

fig.set_tight_layout(True)
fig.savefig("plag_V_KT.pdf")

Expand Down Expand Up @@ -189,7 +213,7 @@
for i in range(7):
ax[i].legend()
ax[i].set_xlabel("$p_{an}$")
ax[i].set_ylabel("modulus (GPa)")
ax[i].set_ylabel("$C_{Nij}$ (GPa)")

fig.set_tight_layout(True)
fig.savefig("plag_stiffnesses.pdf")
Expand All @@ -198,7 +222,7 @@
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$"
name = f"$S_{{N{i+1}{j+1}}} / \\beta_{{NR}}$"
axi = int((irow - (irow) % 3) / 3)
ln = ax[axi].plot(p_ans[mask2], prps["psi"][mask2, i, j])
ax[axi].plot(
Expand All @@ -220,7 +244,7 @@
for i in range(7):
ax[i].legend()
ax[i].set_xlabel("$p_{an}$")
ax[i].set_ylabel("$d\\Psi / df$")
ax[i].set_ylabel("$S_{Nij} / \\beta_{NR}$")

fig.set_tight_layout(True)
fig.savefig("plag_psi.pdf")
Expand Down Expand Up @@ -329,7 +353,7 @@
alpha=0.5,
)

ax[axi].set_ylabel(f"$\\beta_{{{axi}}}$ (m)")
ax[axi].set_ylabel(f"$\\beta_{{T{axi+1}}}$ (GPa$^{{-1}}$)")
ax[axi].set_xlabel("$p_{an}$")

fig.set_tight_layout(True)
Expand Down
Loading