''' ======== Tutorial ======== This notebook demonstrates the basic capabilities of Peridoxia. ***To run this code interactively in your web browser, click the "Launch Binder" button below:*** (Note that the file may take a long time to load, and outputs will not be saved.) .. image:: https://mybinder.org/badge_logo.svg :target: https://mybinder.org/v2/gl/skbirner%2Fperidoxia/main?urlpath=%2Fdoc%2Ftree%2F.%2Fdoc%2Fsource%2Fauto_examples%2Fplot_tutorial.ipynb ''' # %% # Installing and Importing # ======================== # # The first step towards running the model is to ensure that all necessary packages are installed. # If you are running in an interactive MyBinder instance, this has been done for you. # If you are running the code locally, please ensure you have followed the installation steps: https://peridoxia.readthedocs.io/en/latest/installation.html # # Once you have installed the Peridoxia package, you will need to import its modules into your script: # # *NOTE: To run a cell, click inside the shell and hit Shift+Enter. All cells in a notebook can be run from the top "Run" menu ("Run All Cells")* # # *NOTE: You may get some Deprecation Warnings when you import these modules. That's fine!* # --------------------------------------------------------------------------------------------------------------------- # Package Imports: These modules contain the functions and classes necessary for running the model. # --------------------------------------------------------------------------------------------------------------------- from peridoxia import data_structures as data_struct from peridoxia import input_composition as input_comp from peridoxia import main # --------------------------------------------------------------------------------------------------------------------- # General Imports: It's also useful to have a few common packages imported into the code. # --------------------------------------------------------------------------------------------------------------------- import numpy as np import matplotlib.pyplot as plt import pandas as pd # Suppress deprecation warnings in Binder script import warnings warnings.filterwarnings("ignore", category=DeprecationWarning) # %% # Minimum Working Example # ======================= # # This code block demonstrates the minimum amount of code needed to run a. # Then, the code below breaks down the inputs and parameters, demonstrating how to use the model and adapt it to your own use case. pot_TC = 1350 initial_assemblage = input_comp.get_WH05_starting_comp(bulk_Fe2O3=0.3) initial_assemblage_in_garnet_field = input_comp.garnet_project(initial_assemblage, pot_TC, 4.0, spl_final=0.05, display=False) model_run = main.ModelRun( label = 'Tp = 1350°C', potential_temperature_C = pot_TC, initial_assemblage = initial_assemblage_in_garnet_field, ) model_run.plot() # %% # Breaking Down the Code # ====================== # # 1. Determining Starting Composition # ----------------------------------- # # 1a. Using Workman and Hart (2005) starting composition. # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # # Currently, the only 'built-in' input composition is that of DMM from Workman and Hart (2005). # You can access this composition using the following code (change the bulk_Fe2O3 value to your preferred value, if necessary): initial_assemblage = input_comp.get_WH05_starting_comp(bulk_Fe2O3=0.3, display=False) # %% # 1b. Creating your own input composition. # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # # If you wish to use a different starting composition, you can! It just takes a few more steps, # but it's also useful for understanding what the code is doing. # # *Example of how to format an input composition. Change these values as needed.* phase_modes = {'olv':57, 'opx':28, 'cpx':13, 'spl':2, 'gt':0, 'liq':0, 'agg_liq':0} phase_oxides = { # Minerals present in starting composition 'olv': {'SiO2':40.7, 'Al2O3':0, 'Fe2O3':0, 'FeO':10.16, 'MgO':48.59, 'CaO':0.05, 'Cr2O3':0}, 'opx': {'SiO2':53.36, 'Al2O3':6.46, 'Fe2O3':0, 'FeO':6.27, 'MgO':30.55, 'CaO':2.18, 'Cr2O3':0.76}, 'cpx': {'SiO2':50.61, 'Al2O3':7.87, 'Fe2O3':0, 'FeO':2.94, 'MgO':16.19, 'CaO':19.52, 'Cr2O3':1.20}, 'spl': {'SiO2':0, 'Al2O3':57.54, 'Fe2O3':0, 'FeO':12.56, 'MgO':19.27, 'CaO':0, 'Cr2O3':10.23}, # Minerals not present in starting composition 'gt': {'SiO2':0, 'Al2O3':0, 'Fe2O3':0, 'FeO':0, 'MgO':0, 'CaO':0, 'Cr2O3':0}, 'liq': {'SiO2':0, 'Al2O3':0, 'Fe2O3':0, 'FeO':0, 'MgO':0, 'CaO':0, 'Cr2O3':0}, 'agg_liq': {'SiO2':0, 'Al2O3':0, 'Fe2O3':0, 'FeO':0, 'MgO':0, 'CaO':0, 'Cr2O3':0} } # %% # Once you have your phase modes and phase oxides entered above, we will convert it into a PhaseDictionary object, # which will allow us to use some useful methods to adjust the composition for input into the model. phase_dictionaries = data_struct.PhaseDictionaries(phase_oxides, phase_modes) # %% # Next, you can adjust the bulk Fe2O3, if necessary. The method for adjusting bulk Fe2O3 assumes that the FeO reported # in the input above is FeO_total-- both FeO and Fe2O3 -- as might be measured by EPMA. The code adjusts Fe2O3 contents # until the desired bulk Fe2O3 is reached. An important note is that this code only partitions Fe3+ into phases in an # approximate manner. These values will later be re-distributed according to partitioning relationships-- # this is just a first pass. The ratio of Fe3+/tFe values between phases is given by the 'ferric_phases' argument to # the function. By default, spinel will have an Fe3+/tFe ratio twice that of the two pyroxenes. This generally is a good # enough guess for the re-distribution functions to find the actual Fe3+ contents later in the process. phase_dictionaries.adjust_bulk_Fe2O3(Fe2O3_value=0.3, ferric_phases={'opx':1, 'cpx':1, 'spl':2}) # %% # Once you have the composition you like, you will turn it into an Assemblage object. The difference between a PhaseDictionary # and an Assemblage is that an Assemblage object imposes stoichiometric constraints on each phase. This means that oxide # compositions may change slightly from the initial input composition (including the Fe2O3 value indicated above-- I'm currently # writing a function to re-adjust the Fe2O3 back to the desired value.) initial_assemblage = phase_dictionaries.create_assemblage() # This code prints the information for your Assemblage. Make sure it looks right before continuing! print('Assemblage:') initial_assemblage.display_output() print(f"System Fe2O3: {initial_assemblage.get_phase_object('sys').oxides[2]:.5f}") print('NOTE: System Fe2O3 may be slightly different than indicated value, due to adjustments when creating stoichiometric phase objects.') # %% # 2. Projecting Into the Garnet Field # ----------------------------------- # # The next step is to project your input composition from the spinel field into the garnet field, if necessary. This is done by running # subsolidus garnet-out reactions in reverse until spinel is decreased to a certain value (0.05 wt% by default). # The initial assemblage is dependent on the pressure and temperature conditions of your run, so note those below: pot_TC = 1350 PGPa_0 = 4.0 initial_assemblage_in_garnet_field = input_comp.garnet_project(initial_assemblage, pot_TC, PGPa_0, spl_final=0.05, display=False) # You can use the same lines of code as above to check that your garnet-field composition looks reasonable: print('Assemblage in Garnet Field:') initial_assemblage_in_garnet_field.display_output() # %% # 3. Setting Other Input Parameters # --------------------------------- # # The model has other input parameters, most of which can be kept as the defaults, but can be changed if desired. P_step = -0.01 # decrease in pressure, in GPa, at each model step melting = True # toggles melting on and off label = 'Tutorial Run at 1350°C' # label for model run color = 'blue' # color for plotting model run debug = 0 # the level of output to display while the model is running (0, 1,or 2) # %% # 4. Running the Model # -------------------- # # Now you can run the model! The model arguments draw from the information you indicated above, so you don't # need to change them here. Each model run will take about 5 seconds. model_run = main.ModelRun( label = f'Tp = {pot_TC}°C', # this is an "f string" that will insert the value of the variable within the {} brackets into the final string color = color, melting = melting, potential_temperature_C = pot_TC, P_range_GPa = [PGPa_0, 0.1], # ends automatically at cpx-out if melting is enabled P_step = P_step, initial_assemblage = initial_assemblage_in_garnet_field, debug = debug, # change this to 1 or 2 to see detailed output ) # %% # 5. Investigating the Results # ---------------------------- # # 5a. Accessing Output Variables # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # # Model results can be investigated by accessing attributes of the ModelRun object. # The lines below indicate the variables that can be accessed. More info can be found on the API page: https://peridoxia.readthedocs.io/en/latest/api/index.html. # %% # Properties that can be called for the system: # print(model_run.input_arguments) # print(model_run.TK) # print(model_run.TC) # print(model_run.TK_olv_spl) # print(model_run.PGPa) # print(model_run.Pbar) # print(model_run.G_cont_to_fO2) # print(model_run.logfO2_abs) # print(model_run.FFM) # print(model_run.QFM) # print(model_run.logfO2_dFFM) # print(model_run.logfO2_dQFM) # %% # Properties that can be called for each phase: # # *(Here, 'sys' refers to the system, which consists of all solid phases and the incremental liquid, but NOT the aggregated liquid. # To call properties for other phases, use model_run.[phs].[property], replacing [phs] with the phase abbreviation: olv, opx, cpx, spl, gt, liq, agg_liq, # and replacing [property] with the property name.)* # print(model_run.sys.phs_str) # print(model_run.sys.name) # print(model_run.sys.oxide_order) # print(model_run.sys.ext_end_mol) # print(model_run.sys.mass) # print(model_run.sys.oxides) # print(model_run.sys.mol_cat) # print(model_run.sys.Fe3_Al) # print(model_run.sys.Fe3_tFe) # print(model_run.sys.Cr_num) # print(model_run.sys.Mg_num) # print(model_run.olv.Fo_num) # print(model_run.olv.fO2_cont) # print(model_run.opx.Ts_per_6_O) # print(model_run.opx.Al_per_6_O) # print(model_run.opx.XM1XM2) # print(model_run.opx.fO2_cont) # print(model_run.cpx.Ts_per_6_O) # print(model_run.cpx.Al_per_6_O) # print(model_run.spl.aFe3O4) # print(model_run.spl.fO2_cont) # print(model_run.gt.Ca_mass_ratio) # print(model_run.liq.fO2_Kress_and_Carmichael_1991()) # print(model_run.liq.fO2_Borisov_2018(pressure_term='Zhang2017')) # print(model_run.liq.fO2_Hirschmann_2022(pressure_term='Deng2020')) # %% # If you would like a nicely formatted version of a single variable, you can create a dataframe: property = model_run.olv.Fo_num # put the property here # The user doesn't need to change the lines below. df = pd.DataFrame(property) # this line creates a dataframe df_no_idx = df.to_string(index=False,header=False) # this line removes extraneous info # print(df_no_idx) # uncomment this line to display the data! # %% # 5b. Making Plots # ^^^^^^^^^^^^^^^^ # # The easiest way to make a plot of a single run is to use the .plot() method. # By default, the .plot() method creates a plot of logfO2_dFFM vs Pressure. # To make other plots, indicate the x and y values as strings-- any attribute of the ModelRun object is a valid input. model_run.plot() # %% x= 'PGPa' y= 'gt.mass' model_run.plot(x,y) # %% # 5c. Making Data Tables # ^^^^^^^^^^^^^^^^^^^^^^ # # The easiest way to make a table of the main properties of a single run is to use the .get_data_table() method. # You can save the results as a .xlsx by using .save_data_table() method. # Be sure to download your output file to your local computer before exiting the Binder instance! Tables will not be saved! # I also recommend changing the file name below if you have modified the tutorial code. table = model_run.get_data_table() # table # (uncomment to display the table) # %% # To add a column to the table, use the following notation table['logfO2, dQFM'] = model_run.logfO2_dQFM table['Fo#'] = model_run.olv.Fo_num # table # (uncomment to display the table) # %% # To save the table as an Excel file, use the .save_data_table() method: model_run.save_data_table('tutorial_output_table.xlsx') # %% # If you have edited the data table (e.g., added columns), you'll need to specify the new table for saving: model_run.save_data_table('tutorial_output_table_edited.xlsx', data_table=table)