Note
Go to the end to download the full example code.
MagmaForge: Fractional Crystallization#
Demonstrates how to run an equilibrium crystallization routine using MagmaForge.
Open this code in an executable MyBinder instance (MyBinder links may be slow to load– please be patient!):
Initialization#
Import necessary packages:
import thermoengine
from thermoengine import magmaforge
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
Define a bulk composition:
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},)
System and Calculations#
Define system using magmaforge:
sys = magmaforge.System(comp=morb_oxides,
T_liquidus=True,
P_bar=1000.0,
logfO2 = ('QFM', -1),
database='MELTS_v1_0',
)
/workspaces/ThermoEngineLite/thermoengine/thermoengine/magmaforge/system.py:251: UserWarning: Warning: Setting an fO2 value will redistribute the FeO and Fe2O3 values given.
warnings.warn('Warning: Setting an fO2 value will redistribute the FeO and Fe2O3 values given.')
Crystallize system:
sys.crystallize(method='frac',
T_step=4,
fix_fO2=True,
T_final_C = 1000
)
<thermoengine.magmaforge.system.System object at 0xffff8ef5add0>
Output#
Magma Evolution plot
magmaforge.plot.magma_evolution(sys.history)
If you want to specify the color for each phase, you can assign them via a dictionary:
phase_colors = {'Spinel': 'gray',
'Feldspar': 'tan',
'Olivine': 'springgreen',
'Clinopyroxene' : 'green',
'Orthopyroxene': 'darkolivegreen'}
magmaforge.plot.magma_evolution(sys.history, phase_colors=phase_colors)
Return useful output variables#
The following commands demonstrate examples of how to return properties for all states in system’s history. Other properties can be accessed as well– see https://thermoenginelite.readthedocs.io/en/latest/api/magmaforge_system.html for more information.
T = sys.history.get_temperatures(unit='C')
P = sys.history.get_pressures(unit='bar')
logfO2 = sys.history.get_logfO2s(buffer='QFM')
liq_comp = sys.history.get_phase_comps(phase_name='Liquid', unit='wt_oxides')
phase_frac_table = sys.history.get_phase_frac_table(phases='present')
print(phase_frac_table) # replace with whichever variable you would like to access
Olivine Feldspar Clinopyroxene Spinel Liquid
1503.086379 0.000000 0.000000 0.000000 0.000000 1.000000
1499.086379 0.000000 0.000000 0.000000 0.000113 0.999889
1495.086379 0.000000 0.006817 0.000000 0.000093 0.992983
1491.086379 0.002494 0.020345 0.000000 0.000030 0.970181
1487.086379 0.012486 0.034416 0.000000 0.000000 0.923494
1483.086379 0.011316 0.031003 0.000000 0.000000 0.881632
1479.086379 0.010222 0.027934 0.000000 0.000000 0.843887
1475.086379 0.009284 0.025300 0.000000 0.000000 0.809675
1471.086379 0.007684 0.024215 0.008416 0.000000 0.769700
1467.086379 0.005437 0.023586 0.021976 0.000011 0.719094
1463.086379 0.004753 0.020800 0.019368 0.000083 0.674600
1459.086379 0.004226 0.018493 0.016972 0.000077 0.635282
1455.086379 0.003789 0.016565 0.014974 0.000071 0.600280
1451.086379 0.003421 0.014936 0.013293 0.000066 0.568916
1447.086379 0.003111 0.013548 0.011865 0.000063 0.540644
1443.086379 0.002846 0.012356 0.010642 0.000060 0.515023
1439.086379 0.002619 0.011325 0.009588 0.000058 0.491690
1435.086379 0.002423 0.010426 0.008672 0.000057 0.470345
1431.086379 0.002254 0.009638 0.007872 0.000057 0.450737
1427.086379 0.002106 0.008944 0.007170 0.000057 0.432654
1423.086379 0.001978 0.008329 0.006551 0.000059 0.415917
1419.086379 0.001865 0.007782 0.006003 0.000062 0.400371
1415.086379 0.001766 0.007294 0.005518 0.000068 0.385880
1411.086379 0.001678 0.006857 0.005090 0.000078 0.372321
1407.086379 0.001600 0.006467 0.004714 0.000095 0.359579
1403.086379 0.001529 0.006121 0.004393 0.000124 0.347538
1399.086379 0.001463 0.005818 0.004127 0.000171 0.336078
1395.086379 0.001401 0.005557 0.003914 0.000239 0.325081
1391.086379 0.001340 0.005336 0.003750 0.000328 0.314437
1387.086379 0.001280 0.005154 0.003633 0.000440 0.304037
1383.086379 0.001214 0.005020 0.003578 0.000596 0.293734
1379.086379 0.001118 0.005027 0.003752 0.000944 0.283000
1375.086379 0.000833 0.005735 0.005125 0.002365 0.269063
1371.086379 0.000099 0.000643 0.000581 0.000294 0.267590
1367.086379 0.001227 0.010759 0.009580 0.005318 0.240758
1363.086379 0.000693 0.004998 0.003879 0.002168 0.229304
1359.086379 0.000699 0.004547 0.003231 0.001816 0.219141
1355.086379 0.000708 0.004208 0.002750 0.001558 0.210030
1351.086379 0.000715 0.003918 0.002358 0.001350 0.201791
1347.086379 0.000722 0.003665 0.002033 0.001179 0.194284
1343.086379 0.000729 0.003442 0.001760 0.001038 0.187398
1339.086379 0.000736 0.003245 0.001529 0.000920 0.181045
1335.086379 0.000742 0.003069 0.001332 0.000819 0.175152
1331.086379 0.000748 0.002910 0.001164 0.000734 0.169661
1327.086379 0.000753 0.002766 0.001020 0.000661 0.164521
1323.086379 0.000757 0.002634 0.000896 0.000598 0.159694
1319.086379 0.000759 0.002512 0.000789 0.000543 0.155144
1315.086379 0.000759 0.002399 0.000697 0.000495 0.150844
1311.086379 0.000758 0.002294 0.000618 0.000453 0.146768
1307.086379 0.000755 0.002195 0.000550 0.000416 0.142896
1303.086379 0.000751 0.002102 0.000491 0.000383 0.139211
1299.086379 0.000744 0.002015 0.000441 0.000353 0.135698
1295.086379 0.000736 0.001931 0.000397 0.000327 0.132344
1291.086379 0.000727 0.001853 0.000359 0.000304 0.129137
1287.086379 0.000716 0.001778 0.000327 0.000283 0.126067
1283.086379 0.000704 0.001706 0.000299 0.000264 0.123126
1279.086379 0.000691 0.001638 0.000274 0.000248 0.120306
1275.086379 0.000677 0.001573 0.000253 0.000232 0.117601
1271.086379 0.000661 0.001511 0.000235 0.000218 0.115003
If you wish to save the full state history, you can do so with the .save_full_state_history_table() method. Soon, you will be able to reload this table into MagmaForge to access its properties.
sys.history.save_full_state_history_table(filename='magmaforge_output.csv', P_unit='bar', T_unit='K', fO2_buffer='QFM')
Total running time of the script: (0 minutes 18.692 seconds)