Skip to content

Commit ffcc420

Browse files
bobmyhillbobmyhill
committed
refactor benchmarks
1 parent cca74b7 commit ffcc420

File tree

4 files changed

+56
-63
lines changed

4 files changed

+56
-63
lines changed

‎misc/benchmarks/ab_initio_liquids.py‎

Lines changed: 16 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -56,24 +56,29 @@
5656
ax_E.set_xlabel("Volume (cm^3/mol)")
5757

5858
ax_P.set_ylabel("Pressure (GPa)")
59-
ax_S.set_ylabel("Entropy/nR (J)")
59+
ax_S.set_ylabel("Entropy/nR")
6060
ax_E.set_ylabel("Internal Energy (kJ/mol)")
6161

6262
ax_P.imshow(mpimg.imread(plots[i][0][0]), extent=plots[i][0][1], aspect="auto")
6363
ax_S.imshow(mpimg.imread(plots[i][1][0]), extent=plots[i][1][1], aspect="auto")
6464
ax_E.imshow(mpimg.imread(plots[i][2][0]), extent=plots[i][2][1], aspect="auto")
6565

6666
for temperature in temperatures:
67-
pressures = np.empty_like(volumes)
68-
entropies = np.empty_like(volumes)
69-
energies = np.empty_like(volumes)
67+
pressures = np.empty_like(volumes) + np.nan
68+
entropies = np.empty_like(volumes) + np.nan
69+
energies = np.empty_like(volumes) + np.nan
7070

7171
for j, volume in enumerate(volumes):
72-
pressures[j] = phase.method.pressure(temperature, volume, phase.params)
73-
entropies[j] = phase.method._S_el(temperature, volume, phase.params)
74-
energies[j] = phase.method.molar_internal_energy(
75-
0.0, temperature, volume, phase.params
76-
)
72+
try:
73+
phase.set_state_with_volume(volume, temperature)
74+
pressures[j] = phase.pressure
75+
entropies[j] = phase.method._S_el(temperature, volume, phase.params)
76+
energies[j] = phase.method._molar_internal_energy(
77+
phase.pressure, temperature, volume, phase.params
78+
)
79+
energies[j] = phase.molar_internal_energy
80+
except Exception as e:
81+
print(e.message)
7782

7883
ax_P.plot(
7984
volumes * 1.0e6,
@@ -96,6 +101,7 @@
96101

97102
ax_E.legend(loc="upper right")
98103

104+
fig.set_layout_engine("tight")
99105
plt.show()
100106

101107

@@ -144,4 +150,5 @@
144150
ax_hugoniot.set_ylabel("Pressures (GPa)")
145151
ax_hugoniot.set_title("hugoniot")
146152

153+
fig.set_layout_engine("tight")
147154
plt.show()

‎misc/benchmarks/benchmark.py‎

Lines changed: 29 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -453,40 +453,35 @@ def check_slb_fig3():
453453
"""
454454
Benchmark grueneisen parameter against figure 3 of Stixrude and Lithgow-Bertelloni (2005b)
455455
"""
456-
perovskite = burnman.Mineral()
457-
perovskite.params = {
458-
"name": "perovksite",
459-
"V_0": molar_volume_from_unit_cell_volume(168.27, 4.0),
460-
"grueneisen_0": 1.63,
461-
"q_0": 1.7,
462-
}
463-
464-
volume = np.linspace(0.6, 1.0, 100)
465-
grueneisen_slb = np.empty_like(volume)
466-
grueneisen_mgd = np.empty_like(volume)
467-
q_slb = np.empty_like(volume)
468-
q_mgd = np.empty_like(volume)
456+
perovskite = burnman.minerals.SLB_2005.mg_perovskite()
457+
perovskite.params["V_0"] = molar_volume_from_unit_cell_volume(168.27, 4.0)
458+
perovskite.params["grueneisen_0"] = 1.63
459+
perovskite.params["q_0"] = 1.7
460+
perovskite.set_method(slb.SLB2)
461+
462+
Vrel = np.linspace(0.6, 1.0, 100)
463+
grueneisen_slb = np.empty_like(Vrel)
464+
grueneisen_mgd = np.empty_like(Vrel)
465+
q_slb = np.empty_like(Vrel)
466+
q_mgd = np.empty_like(Vrel)
469467

470-
slb_eos = slb.SLB2()
471468
mgd_eos = mgd.MGD2()
472469

473470
# calculate its thermal properties
474-
for i in range(len(volume)):
475-
# call with dummy pressure and temperatures, they do not change it
476-
grueneisen_slb[i] = slb_eos.grueneisen_parameter(
477-
0.0, 0.0, volume[i] * perovskite.params["V_0"], perovskite.params
478-
)
479-
grueneisen_mgd[i] = mgd_eos.grueneisen_parameter(
480-
0.0, 0.0, volume[i] * perovskite.params["V_0"], perovskite.params
481-
)
482-
q_slb[i] = slb_eos.volume_dependent_q(1.0 / volume[i], perovskite.params)
471+
for i in range(len(Vrel)):
472+
volume = Vrel[i] * perovskite.params["V_0"]
473+
perovskite.set_state_with_volume(volume, 300.0)
474+
grueneisen_slb[i] = perovskite.grueneisen_parameter
475+
476+
grueneisen_mgd[i] = mgd_eos._grueneisen_parameter(volume, perovskite.params)
477+
q_slb[i] = slb.SLB2()._volume_dependent_q(1.0 / Vrel[i], perovskite.params)
483478
q_mgd[i] = perovskite.params["q_0"]
484479

485480
# compare with figure 7
486481
fig1 = mpimg.imread("../../burnman/data/input_figures/slb_fig3.png")
487482
plt.imshow(fig1, extent=[0.6, 1.0, 0.35, 2.0], aspect="auto")
488-
plt.plot(volume, grueneisen_slb, "g+", volume, grueneisen_mgd, "b+")
489-
plt.plot(volume, q_slb, "g+", volume, q_mgd, "b+")
483+
plt.plot(Vrel, grueneisen_slb, "g+", Vrel, grueneisen_mgd, "b+")
484+
plt.plot(Vrel, q_slb, "g+", Vrel, q_mgd, "b+")
490485
plt.xlim(0.6, 1.0)
491486
plt.ylim(0.35, 2.0)
492487
plt.ylabel("Grueneisen parameter")
@@ -517,37 +512,28 @@ def check_slb_fig7_txt():
517512
}
518513
forsterite.set_method("slb3")
519514

520-
data = np.loadtxt("../../burnman/data/input_minphys/slb_fig7.txt", skiprows=2)
515+
data = np.loadtxt(
516+
"../../burnman/data/input_minphys/slb_fig7.txt", skiprows=2, unpack=True
517+
)
521518

522-
temperature = np.array(data[:, 2])
523-
pressure = np.array(data[:, 0])
524-
rho = np.array(data[:, 3])
519+
pressure, _, temperature = data[:3]
520+
temperature[0] = 1.0e-1
521+
rho, Kt, Ks, G, VB, VS, VP = data[3:10]
525522
rho_comp = np.empty_like(rho)
526-
Kt = np.array(data[:, 4])
527523
Kt_comp = np.empty_like(Kt)
528-
Ks = np.array(data[:, 5])
529524
Ks_comp = np.empty_like(Ks)
530-
G = np.array(data[:, 6])
531525
G_comp = np.empty_like(G)
532-
VB = np.array(data[:, 7])
533526
VB_comp = np.empty_like(VB)
534-
VS = np.array(data[:, 8])
535527
VS_comp = np.empty_like(VS)
536-
VP = np.array(data[:, 9])
537528
VP_comp = np.empty_like(VP)
538-
vol = np.array(data[:, 10])
529+
530+
vol, alpha, Cp, gr, gibbs, entropy, enthalpy = data[10:17]
539531
vol_comp = np.empty_like(vol)
540-
alpha = np.array(data[:, 11])
541532
alpha_comp = np.empty_like(alpha)
542-
Cp = np.array(data[:, 12])
543533
Cp_comp = np.empty_like(Cp)
544-
gr = np.array(data[:, 13])
545534
gr_comp = np.empty_like(gr)
546-
gibbs = np.array(data[:, 14])
547535
gibbs_comp = np.empty_like(gibbs)
548-
entropy = np.array(data[:, 15])
549-
entropy_comp = np.empty_like(gibbs)
550-
enthalpy = np.array(data[:, 16])
536+
entropy_comp = np.empty_like(entropy)
551537
enthalpy_comp = np.empty_like(gibbs)
552538

553539
for i in range(len(temperature)):

‎misc/benchmarks/liquid_iron_comparison.py‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@
1616
liq, P=10.0e9, T=7000.0, tol=1.0e-3, verbose=True, including_shear_properties=False
1717
)
1818

19+
R = burnman.constants.gas_constant
1920

2021
# Find heat capacities
2122
temperatures = np.linspace(1000.0, 15000.0, 101)
@@ -26,10 +27,8 @@
2627
for rho in densities:
2728
V = m / rho
2829
for i, T in enumerate(temperatures):
29-
Cvs[i] = (
30-
liq.method.molar_heat_capacity_v(0.0, T, V, liq.params)
31-
/ burnman.constants.gas_constant
32-
)
30+
liq.set_state_with_volume(V, T)
31+
Cvs[i] = liq.molar_heat_capacity_v / R
3332

3433
plt.plot(temperatures, Cvs)
3534

@@ -60,6 +59,7 @@ def isentrope(T, P, Sref, mineral):
6059
rhos = np.empty_like(pressures)
6160
Vps = np.empty_like(pressures)
6261

62+
6363
for i, P in enumerate(pressures):
6464
temperatures[i] = fsolve(isentrope, 1811.0, args=(P, reference_entropy, liq))[0]
6565
rhos[i] = liq.density

‎misc/ref/calibrant_benchmarks.py.out‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -31,12 +31,12 @@ V=1.9044e-05 m^3/mol (last row):
3131

3232
Checking burnman.calibrants.Dorogokupets_2007.Ag...
3333
V=0.7*V0, P=100.37 GPa, T=298.15 K: 0.250 MPa difference
34-
V=0.7*V0, P=104.73 GPa, T=1000 K: -0.662 MPa difference
34+
V=0.7*V0, P=104.73 GPa, T=1000 K: -0.661 MPa difference
3535
V=0.7*V0, P=111.15 GPa, T=2000 K: -0.997 MPa difference
3636
Checking burnman.calibrants.Dorogokupets_2007.Al...
3737
V=0.7*V0, P=56.97 GPa, T=298.15 K: -11.383 MPa difference
38-
V=0.7*V0, P=60.03 GPa, T=1000 K: -12.426 MPa difference
39-
V=0.7*V0, P=64.86 GPa, T=2000 K: -43.175 MPa difference
38+
V=0.7*V0, P=60.03 GPa, T=1000 K: -12.427 MPa difference
39+
V=0.7*V0, P=64.86 GPa, T=2000 K: -43.176 MPa difference
4040
Checking burnman.calibrants.Dorogokupets_2007.Au...
4141
V=0.7*V0, P=164.83 GPa, T=298.15 K: -0.482 MPa difference
4242
V=0.7*V0, P=169.00 GPa, T=1000 K: -4.064 MPa difference
@@ -59,11 +59,11 @@ V=0.7*V0, P=224.83 GPa, T=1000 K: -19.063 MPa difference
5959
V=0.7*V0, P=228.53 GPa, T=2000 K: -20.886 MPa difference
6060
Checking burnman.calibrants.Dorogokupets_2007.MgO...
6161
V=0.7*V0, P=116.72 GPa, T=298.15 K: 52.963 MPa difference
62-
V=0.7*V0, P=120.96 GPa, T=1000 K: 51.888 MPa difference
63-
V=0.7*V0, P=128.13 GPa, T=2000 K: 51.510 MPa difference
62+
V=0.7*V0, P=120.96 GPa, T=1000 K: 51.887 MPa difference
63+
V=0.7*V0, P=128.13 GPa, T=2000 K: 51.509 MPa difference
6464
Checking burnman.calibrants.Dorogokupets_2007.diamond...
65-
V=0.7*V0, P=301.53 GPa, T=298.15 K: 13.945 MPa difference
66-
V=0.7*V0, P=303.90 GPa, T=1000 K: 13.360 MPa difference
65+
V=0.7*V0, P=301.53 GPa, T=298.15 K: 13.944 MPa difference
66+
V=0.7*V0, P=303.90 GPa, T=1000 K: 13.359 MPa difference
6767
V=0.7*V0, P=309.77 GPa, T=2000 K: 10.797 MPa difference
6868

6969
Checking burnman.calibrants.Tsuchiya_2003.Au at compression 0.34...

0 commit comments

Comments
 (0)