Hello everyone,
I am currently learning to use alphaMELTS by following the examples provided on the ENKI website. Two specific examples demonstrate modeling systematic variation for the same initial conditions:
Fractionation example: 1. https://enki-portal.gitlab.io/ThermoEngine/MELTS-v1.0.2-fractionation.html
Equilibrium example: 2. https://enki-portal.gitlab.io/ThermoEngine/MELTS-v1.0.2-equilibrium.html
To deepen my understanding, I tried translating these examples into Python scripts. However, I encountered issues running the Python code. Below is one of my attempts:
~~~~~~~~
from thermoengine import equilibrate
import matplotlib.pyplot as plt
import numpy as np
# Initialize MELTS v1.0.2 instance
melts = equilibrate.MELTSmodel()
# Set initial bulk composition (wt% or grams)
feasible = melts.set_bulk_composition({
'SiO2': 77.5, 'TiO2': 0.08, 'Al2O3': 12.5, 'Fe2O3': 0.207,
'Cr2O3': 0.0, 'FeO': 0.473, 'MnO': 0.0, 'MgO': 0.03,
'NiO': 0.0, 'CoO': 0.0, 'CaO': 0.43, 'Na2O': 3.98,
'K2O': 4.88, 'P2O5': 0.0, 'H2O': 5.5
})
# Optional: Suppress unwanted phases
melts.set_phase_inclusion_status({'Nepheline': False, 'OrthoOxide': False})
# Compute equilibrium at specific T (°C) and P (MPa)
output = melts.equilibrate_tp(760.0, 175.0, initialize=True)
(status, t, p, xmlout) = output
print(status, t, p)
# Display summary output
melts.output_summary(xmlout)
# List thermodynamic properties for the system
props = melts.get_list_of_properties()
for prop in props:
print(f"{prop:<20s} {melts.get_property_of_phase(xmlout, 'System', prop):13.6e} {melts.get_units_of_property(prop):<10s}")
# List chemical affinities and undersaturated phase compositions
affinity_dict = melts.get_dictionary_of_affinities(xmlout, sort=True)
for phase in affinity_dict:
affinity, formulae = affinity_dict[phase]
print(f"{phase:<20s} {affinity:10.2f} {formulae:<60s}")
# Set up for sequential calculations and plotting
number_of_steps = 20
t_increment = -1.0
p_increment = 0.0
plotPhases = ['Liquid', 'Sanidine', 'Plagioclase', 'Quartz', 'Water']
plotColors = ['ro', 'bo', 'go', 'co', 'mo']
# Initialize data arrays for plotting
n = len(plotPhases)
xPlot = np.zeros(number_of_steps + 1)
yPlot = np.zeros((n, number_of_steps + 1))
xPlot
for i in range(n):
yPlot
- = melts.get_property_of_phase(xmlout, plotPhases)
# Create a plot for the phase data
plt.ion()
fig, ax = plt.subplots()
ax.set_xlim([min(t, t + t_increment * number_of_steps), max(t, t + t_increment * number_of_steps)])
ax.set_ylim([0.0, 100.0])
graphs = [ax.plot(xPlot, yPlot, plotColors)
ax.legend(graphs, plotPhases, loc='upper left')
# Perform sequential calculations and update the plot
for i in range(1, number_of_steps + 1):
output = melts.equilibrate_tp(t + i * t_increment, p + i * p_increment)
status, t, p, xmlout = output
print(f"{status:<30s} {t:8.2f} {p:8.2f}")
xPlot = t
for j in range(n):
yPlot[j] = melts.get_property_of_phase(xmlout, plotPhases[j])
for j, graph in enumerate(graphs):
graph.set_xdata(xPlot)
graph.set_ydata(yPlot[j])
fig.canvas.draw()
# Show the plot (important for standard Python scripts)
plt.show()
~~~~~~~~~~~
I know these examples are designed to run on the ENKI server. However, I encountered difficulties translating them into a Python script, particularly because I couldn't find the thermoengine module referenced in the original examples. This has been a blocking issue, preventing me from running the code successfully on my local machine.
Could someone kindly help me verify the code or guide me in correctly setting it up and running it in Python, so I can run it smoothly, just like the example from this web link (https://github.com/magmasource/alphaMELTS/blob/main/examples/py/tutorial.py), using the "alphamelts-py-2.3.1-win64" package on my personal computer?
I would greatly appreciate any advice or insights.
Thank you in advance for your help!
Best regards,
Shanke
It should be ok now.
Paula