MagmaForge: PhaseMap

Demonstrates how to find phase assemblages over a grid of P/T conditions.

Open this code in an executable MyBinder instance (MyBinder links may be slow to load– please be patient!):

Initialization

Import necessary packages:

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},)

Run crystallization paths at a set of pressures

P_kbar = [10,9,8,7,6,5,4,3,2,1] # to get highest T necessary first
phase_frac_tables = []

max_len = 0
for idx,P in enumerate(P_kbar):

    if idx==0:
        sys = magmaforge.System(comp=morb_oxides,
                            T_liquidus = True,
                            P_kbar=P,
                            database='MELTS_v1_0',
                            )

        T_max = sys.T_C

    else:
        sys = magmaforge.System(comp=morb_oxides,
                            T_C = T_max,
                            P_kbar=P,
                            database='MELTS_v1_0',
                            )

    sys.crystallize(method='equil',
                    T_step=100,
                    fix_fO2=False,
                    T_final_C = 600
                    )

    tbl = sys.history.get_phase_frac_table('all')

    if len(tbl.index) > max_len:
        max_len = len(tbl.index)

    phase_frac_tables.append(tbl)

Plot phase fractions

fig, axs = plt.subplots(len(P_kbar), max_len, figsize=(6,6))

# Olivine  Feldspar  Clinopyroxene  Orthopyroxene    Spinel Quartz  Fayalite  Magnetite  Tridymite   Kyanite  Andalusite Sillimanite  Water  Liquid
color_series = pd.Series('k', index=phase_frac_tables[0].columns)
color_series['Liquid'] = 'red'
color_series['Water'] = 'blue'
color_series['Spinel'] = 'gray'
color_series['Feldspar'] = 'tan'
color_series['Olivine'] = 'springgreen'
color_series['Clinopyroxene'] = 'green'
color_series['Clinopyroxene_1'] = 'green'
color_series['Orthopyroxene'] = 'darkolivegreen'

for i, table in enumerate(phase_frac_tables):
    imax = len(phase_frac_tables)-1
    for j, (index, row) in enumerate(table.iterrows()):
        #jmax = len(phase_frac_tables)-1
        if i==0:
            axs[imax,j].set_xlabel(f'{index-273.15:.0f}')
        color_list = []
        for phase in row.index:
            if phase in color_series.index:
                color_list.append(color_series[phase])
            else:
                color_list.append('pink')
        axs[imax-i,j].pie(row, colors=color_list)
        axs[imax-i,j].axis('equal')

    while j+1<max_len:
        j+=1
        axs[imax-i,j].pie([1], colors=['white'])

    axs[i,0].set_ylabel(P_kbar[imax-i])

fig.supxlabel("Temperature (°C)")
fig.supylabel("Pressure (kbar)")
plt.tight_layout()

plt.show()