Note
Go to the end to download the full example code.
Tutorial: MagmaForge#
Demonstrates how to use MagmaForge to perform MELTS-style calculations.
Open this code in an executable MyBinder instance (MyBinder links may be slow to load– please be patient!):
Introduction#
MagmaForge is a tool for easily running MELTS-style calculations such as progressive crystallization of liquids. It is designed to be a user-friendly wrapper for the underlying phase model code and equillibration routines available in ThermoEngine. This tutorial demonstrates some of the capabilities available in MagmaForge; it is a work in progress and will be updated over the next few months to demonstrate more capabilities.
Minimum Working Example#
The following code block demonstrates the minimum amount of code necessary to run an equilibrium crystallization routine using MagmaForge. In the following sections, we will go through each part and demonstrate additional capabilities.
from thermoengine import magmaforge
import pandas as pd # a useful data analysis package
morb_oxides = pd.Series({
'SiO2': 48.68,
'TiO2': 1.01,
'Al2O3': 17.64,
'Fe2O3': 0.89,
'Cr2O3': 0.0425,
'FeO': 7.59,
'MgO': 9.10,
'CaO': 12.45,
'Na2O': 2.65,
'K2O': 0.03,
'P2O5': 0.08,
'H2O': 0.2},) # Values in grams (extensive units)
sys = magmaforge.System(comp = morb_oxides,
P_bar = 1000.0, # bars
T_K = 1600, # K
logfO2 = ('NNO',-1), # setting fO2 equal to 1 log unit below the NNO buffer
database='MELTS_v1_0', # liquid model name
)
sys.crystallize(
method='equil', # equilibrium crystallization
Tstep=5, # decrease in temperature at each step, in K
fix_fO2=True, # hold fO2 constant at initial fO2 value, relative to NNO (more buffers to be implemented soon!)
)
magmaforge.plot.phase_fractions(sys.history)
/workspaces/ThermoEngineLite/thermoengine/thermoengine/magmaforge/system.py:265: UserWarning: Warning! Setting an fO2 value will overwrite the FeO and Fe2O3 values given.
warnings.warn('Warning! Setting an fO2 value will overwrite the FeO and Fe2O3 values given.')
Breaking It Down#
1. Initialization#
To use MagmaForge, you’ll first need to import the ThermoEngine package and the MagmaForge module. If you’re running the code locally, you should have already installed ThermoEngine to your computer– here you’re importing it from your computer’s package library into this notebook. If you’re running the code via MyBinder, ThermoEngine is already installed by default in the virtual environment.
from thermoengine import magmaforge
Import Pandas (https://pandas.pydata.org/), a package which provides convenient capabilities for entering, storing, and anaylzing data. If you used Anaconda to create a virtual environment on your machine, Pandas should exist on your machine and the import statement below should work. If you get an error message, you can install Pandas by typing “pip install pandas” into your command line.
import pandas as pd
2. Defining a System#
The first step for most MELTS-style calculations is to define a chemical system– i.e., a bulk chemical composition equilibrated at a set of conditions (usually pressure and temperature, but sometimes other variables like entropy or volume). Here, let’s imagine that we’re interested in a mid-ocean ridge basalt. We can define the bulk chemical composition using a series, which is a data type specific to the Pandas package– the ‘pd’ indicates that we are invoking this package.
morb_oxides = pd.Series({
'SiO2': 48.68,
'TiO2': 1.01,
'Al2O3': 17.64,
'Fe2O3': 0.89,
'Cr2O3': 0.0425,
'FeO': 7.59,
'MgO': 9.10,
'CaO': 12.45,
'Na2O': 2.65,
'K2O': 0.03,
'P2O5': 0.08,
'H2O': 0.2},) # values in grams (extensive units)
Now that we have a bulk chemical composition, we can use MagmaForge to equilibrate the MORB composition at the conditions of interest. We use the magmaforge.System() call to create a system object at the indicated composition, pressure, temperature, and fO2. Here, the initial temperature is chosen to be above the liquidus of the system so that later calculations can crystallize the system.
sys = magmaforge.System(comp = morb_oxides,
P_bar = 1000.0, # bars
T_K = 1600, # K
logfO2 = ('NNO', -1), # setting fO2 equal to 1 log unit below the NNO buffer
database='MELTS_v1_0', # liquid model name
)
/workspaces/ThermoEngineLite/thermoengine/thermoengine/magmaforge/system.py:265: UserWarning: Warning! Setting an fO2 value will overwrite the FeO and Fe2O3 values given.
warnings.warn('Warning! Setting an fO2 value will overwrite the FeO and Fe2O3 values given.')
The sys object that we created in the last cell represents our equilibrated chemical system. We can probe this system to ensure our calculation had the expected behavior. For instance, in the cell below, we are querying the temperature (in K), pressure (in bars), and melt fraction of the system.
print(sys.T_C)
print(sys.P_GPa)
print(sys.melt_fraction)
1326.85
0.1
1.0
3. Performing a Calculation Along a Path#
Now that we have our chemical system, we can perform calculations such as crystallization during cooling. To do this, we will apply the crystallize method to our system. We tell it that we want to use the ‘equil’ (equilibrium/batch crystallization) method of crystallizing, that we want to decrease in temperature steps of 15°, and that we want to hold fO2 steady at the value we set for the system– here NNO -1.
This calculation may return a message reading “No check made for phase separation.” This is fine!
sys.crystallize(
method='equil', # equilibrium crystallization
Tstep=15, # decrease in temperature at each step, in K
fix_fO2=True, # hold fO2 constant at initial fO2 value, relative to NNO (more buffers to be implemented soon!)
)
Bad phase composition for Liquid: [ 9.51536524e-08 3.16490592e-06 -7.97877731e-07 4.12861483e-07
3.08050640e-09 -4.63763903e-04 0.00000000e+00 6.90490869e-05
0.00000000e+00 0.00000000e+00 5.63087312e-07 0.00000000e+00
1.59575585e-06 5.61144607e-04 3.83830588e-05]
/workspaces/ThermoEngineLite/thermoengine/thermoengine/magmaforge/system.py:618: UserWarning: Warning: Unsuccessful cooling step ended equilibrium calculation.
warnings.warn("Warning: Unsuccessful cooling step ended equilibrium calculation.")
<thermoengine.magmaforge.system.System object at 0xffff936f36a0>
Accessing/Plotting Output#
1. Accessing Output as a Plot#
Once you have performed your calculations, you probably want to see the results. MagmaForge has a few built-in plotting functions for creating basic plots. MagmaForge has a few simple built-in plots that can be called using the magmaforge.plot submodule. Here we’re using the .phase_fractions() method to create an automatic plot of phase fractions with respect to temperature as the system cools. The sys.history object contains the history of the system’s state after each calculation and will be very useful for accessing output.
magmaforge.plot.phase_fractions(sys.history)
For a somewhat more in-depth look at our results, we can use magmaforge.plot.magma_evolution() to plot the phase fractions and magma composition during crystallization.
magmaforge.plot.magma_evolution(sys.history)
You can, of course, also create custom plots using the output of the run. We’ll show some examples in a later section.
2. Accessing Output as a Table#
You may also want to access numerical output for your run. Recall that the sys.history object stores the information of the state of the system after each calculation. We can access a table showing the fractions of each phase after each calculation step using the .get_phase_frac_table() method. Here we’re using the ‘present’ argument to only display columns for phases that are present at some point during the run.
sys.history.get_phase_frac_table(phases='present')
The composition of the liquid can be accessed using the .get_liquid_comp_table() method:
sys.history.get_liquid_comp_table()
The composition of phases can be accessed using the .get_mineral_comp_table() method:
sys.history.get_mineral_comp_table('Feldspar')
Any of these tables can be easily downloaded as a csv file using built-in Pandas commands:
If you are running in a MyBinder, be sure to download your output file to your machine before closing your session!
table_to_export = sys.history.get_phase_frac_table(phases='present')
table_to_export.to_csv('output.csv', index=True)
3. Accessing outputs as variables/data structures#
You can also access each property of the system as an array or other data structure, which can be helpful if you want to do further math with the outputs.
T = sys.history.get_temperatures('K')
P = sys.history.get_pressures('bar')
S = sys.history.get_total_entropies()
F = sys.history.get_instantaneous_melt_fractions()
phase_frac = sys.history.get_phase_fractions()
comp_liq = sys.history.get_liquid_comps()
comp_phs = sys.history.get_mineral_comps('Feldspar')
print(T) # replace with whichever variable you'd like to view
[1600. 1585. 1570. 1555. 1540. 1525. 1510. 1495. 1480. 1465. 1450. 1435.
1420. 1405. 1390. 1375. 1360. 1345. 1330. 1315. 1300. 1285. 1270. 1255.
1240.]
More methods are documented in the MagmaForge API reference: https://thermoenginelite.readthedocs.io/en/latest/api/index.html
Total running time of the script: (0 minutes 23.765 seconds)