Source code for VESIcal.models.magmasat

from VESIcal import calibration_checks
from VESIcal import core
from VESIcal import model_classes
from VESIcal import sample_class
from VESIcal import vplot
from VESIcal import batchfile  # needed for status_bar functions

from copy import deepcopy
import numpy as np
import pandas as pd
import warnings as w
import sys

from contextlib import redirect_stdout
import io

# Variable to send MagmaSat warnings into the void
_f = io.StringIO()

# optional dependency thermoengine
try:
    from thermoengine import equilibrate
except ImportError:
    w.warn("\n"
           "\n WARNING: Thermoengine is not installed. MagmaSat model will not function. A "
           "model name must be given when performing any calculation as "
           "model=\'some-model-name\'.", stacklevel=2)
    equilibrate = None

# filter warnings
w.filterwarnings("ignore", category=DeprecationWarning)
w.filterwarnings("ignore", message="rubicon.objc.ctypes_patch has only been tested ")
w.filterwarnings("ignore", message="The handle")


[docs] class MagmaSat(model_classes.Model): """ An object to instantiate a thermoengine equilibrate class """ def __init__(self): if equilibrate is None: raise core.Error( "\n" "ERROR: Thermoengine is not installed. MagmaSat model will not function. A " "model name must be given when performing any calculation as " "model=\'some-model-name\'." ) # -------------- MELTS preamble --------------- # # instantiate thermoengine equilibrate MELTS instance melts = equilibrate.MELTSmodel("1.2.0") # Suppress phases not required in the melts simulation phases = melts.get_phase_names() for phase in phases: melts.set_phase_inclusion_status({phase: False}) melts.set_phase_inclusion_status({"Fluid": True, "Liquid": True}) self.melts = melts # --------------------------------------------- # self.melts_version = ( "1.2.0" # just here so users can see which version is being used ) self.set_volatile_species(["H2O", "CO2"]) self.set_calibration_ranges( [ calibration_checks.CalibrationRange( "pressure", [0.0, 20000.0], calibration_checks.crf_Between, "bar", "MagmaSat", fail_msg=calibration_checks.crmsg_Between_fail, pass_msg=calibration_checks.crmsg_Between_pass, description_msg=calibration_checks.crmsg_Between_description, ), calibration_checks.CalibrationRange( "temperature", [800, 1400], calibration_checks.crf_Between, "oC", "MagmaSat", fail_msg=calibration_checks.crmsg_Between_fail, pass_msg=calibration_checks.crmsg_Between_pass, description_msg=calibration_checks.crmsg_Between_description, ), ] ) self.model_type = "MagmaSat"
[docs] def preprocess_sample(self, sample): """ Returns sample with 0.0 values for any oxides not passed. Parameters ---------- sample: Sample class Magma major element composition. Returns ------- Sample class object """ _sample = deepcopy(sample) """ A "fixedvolatiles" normalization must be applied here to achieve consistent behavior with MagmaSat and between VESIcal calculations. See pull reqeust from 23 June 2022 for discussion of this issue. """ _sample = _sample.get_composition( units="wtpt_oxides", normalization="fixedvolatiles", asSampleClass=True, ) # set any necessary magmasat oxides to 0 if no value given for oxide in core.magmasat_oxides: if oxide in _sample.get_composition(): pass else: _sample.change_composition({oxide: 0.0}) # remove any non-magmasat oxides for oxide in _sample.get_composition().index: if oxide in core.magmasat_oxides: pass else: _sample.delete_oxide(oxide) self.bulk_comp_orig = _sample.get_composition(units="wtpt_oxides") return _sample
[docs] def check_calibration_range(self, parameters, **kwargs): """Checks whether supplied parameters and calculated results are within the calibration range of the model, defined by the CalibrationRange objects. An empty string will be returned if all parameters are within the calibration range. If a parameter is not within the calibration range, a description of the problem will be returned in the string. Parameters ---------- parameters dict Dictionary keys are the names of the parameters to be checked, e.g., pressure temperature, SiO2, etc. Values are the values of each parameter. A complete set need not be given. Returns ------- str String description of any parameters falling outside of the calibration range. """ s = "" for cr in self.calibration_ranges: if cr.check(parameters) is False: s += cr.string(parameters, report_nonexistance=False) if "notsaturated" in kwargs: s += "Sample not saturated at these conditions." return s
[docs] def get_calibration_range(self): """Returns a string describing the calibration ranges defined by the CalibrationRange objects for the model. Returns ------- str String description of the calibration range objects.""" s = "" for cr in self.calibration_ranges: s += cr.string(None) return s
[docs] def get_fluid_mass(self, sample, temperature, pressure, H2O, CO2): """An internally used function to calculate fluid mass. Parameters ---------- sample: Sample class Magma major element composition. temperature: float Temperature in degrees C. pressure: float Pressure in bars H2O: float wt% H2O in the system CO2: float wt% CO2 in the system Returns ------- float mass of the fluid in grams """ pressureMPa = pressure / 10.0 bulk_comp_dict = sample.get_composition(units="wtpt_oxides") bulk_comp_dict = {oxide: bulk_comp_dict[oxide] for oxide in core.magmasat_oxides} bulk_comp_dict["H2O"] = H2O bulk_comp_dict["CO2"] = CO2 self.melts.set_bulk_composition(bulk_comp_dict) with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, pressureMPa, initialize=True) (status, temperature, pressureMPa, xmlout) = output[0] fluid_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Fluid") return fluid_mass
[docs] def get_XH2O_fluid(self, sample, temperature, pressure, H2O, CO2): """An internally used function to calculate fluid composition. Parameters ---------- sample: Sample class Magma major element composition. temperature: float Temperature in degrees C. pressure: float Pressure in bars H2O: float wt% H2O in the system CO2: float wt% CO2 in the system Returns ------- float Mole fraction of H2O in the H2O-CO2 fluid """ _sample = self.preprocess_sample(sample) pressureMPa = pressure / 10.0 _sample.change_composition({"H2O": H2O, "CO2": CO2}) self.melts.set_bulk_composition( _sample.get_composition(units="wtpt_oxides", normalization="none") ) with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, pressureMPa, initialize=True) (status, temperature, pressureMPa, xmlout) = output[0] fluid_comp = self.melts.get_composition_of_phase( xmlout, phase_name="Fluid", mode="component" ) # NOTE mode='component' returns endmember component keys with values in mol fraction. if "Water" in fluid_comp: H2O_fl = fluid_comp["Water"] else: H2O_fl = 0.0 return H2O_fl
[docs] def calculate_dissolved_volatiles( self, sample, temperature, pressure, X_fluid=1, H2O_guess=0.0, verbose=False, **kwargs ): """ Calculates the amount of H2O and CO2 dissolved in a magma at saturation at the given P/T conditions and fluid composition. Fluid composition will be matched to within 0.0001 mole fraction. Parameters ---------- sample: Sample class Magma major element composition. temperature: float or int Temperature, in degrees C. presure: float or int Pressure, in bars. X_fluid: float or int The default value is 1. The mole fraction of H2O in the H2O-CO2 fluid. X_fluid=1 is a pure H2O fluid. X_fluid=0 is a pure CO2 fluid. verbose: bool OPTIONAL: Default is False. If set to True, returns H2O and CO2 concentration in the melt, H2O and CO2 concentration in the fluid, mass of the fluid in grams, and proportion of fluid in the system in wt%. Returns ------- dict A dictionary of dissolved volatile concentrations in wt% with keys H2O and CO2. """ _sample = self.preprocess_sample(sample) if isinstance(X_fluid, int) or isinstance(X_fluid, float): pass else: raise core.InputError("X_fluid must be type int or float") if isinstance(H2O_guess, int) or isinstance(H2O_guess, float): pass else: raise core.InputError("H2O_guess must be type int or float") pressureMPa = pressure / 10.0 if X_fluid != 0 and X_fluid != 1: if X_fluid < 0.001 or X_fluid > 0.999: raise core.InputError( "X_fluid is calculated to a precision of 0.0001 mole " "fraction. Value for X_fluid must be between 0.0001 and " "0.9999." ) H2O_val = H2O_guess CO2_val = 0.0 fluid_mass = 0.0 while fluid_mass <= 0: if X_fluid == 0: CO2_val += 0.1 elif X_fluid >= 0.5: H2O_val += 0.2 # NOTE this is setting XH2Owt of the system (not of the fluid) to X_fluid CO2_val = (H2O_val / X_fluid) - H2O_val # TODO this is what needs to be higher for higher XH2O. Slows down computation # by a second or two else: H2O_val += 0.1 # NOTE this is setting XH2Owt of the system (not of the fluid) to X_fluid CO2_val = (H2O_val / X_fluid) - H2O_val # TODO this is what needs to be higher for higher XH2O. Slows down computation # by a second or two fluid_mass = self.get_fluid_mass( _sample, temperature, pressure, H2O_val, CO2_val ) _sample.change_composition( {"H2O": H2O_val, "CO2": CO2_val}, units="wtpt_oxides" ) self.melts.set_bulk_composition( _sample.get_composition(units="wtpt_oxides", normalization="none") ) with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, pressureMPa, initialize=True) (status, temperature, pressureMPa, xmlout) = output[0] liquid_comp = self.melts.get_composition_of_phase( xmlout, phase_name="Liquid", mode="oxide_wt" ) fluid_comp = self.melts.get_composition_of_phase( xmlout, phase_name="Fluid", mode="component" ) if "Water" in fluid_comp: H2O_fl = fluid_comp["Water"] else: H2O_fl = 0.0 XH2O_fluid = H2O_fl # ------ Coarse Check ------ # while XH2O_fluid < X_fluid - 0.1: # too low coarse check H2O_val += 0.2 XH2O_fluid = self.get_XH2O_fluid( sample, temperature, pressure, H2O_val, CO2_val ) while XH2O_fluid > X_fluid + 0.1: # too high coarse check CO2_val += 0.1 XH2O_fluid = self.get_XH2O_fluid( sample, temperature, pressure, H2O_val, CO2_val ) # ------ Refinement 1 ------ # while XH2O_fluid < X_fluid - 0.01: # too low refinement 1 H2O_val += 0.05 XH2O_fluid = self.get_XH2O_fluid( sample, temperature, pressure, H2O_val, CO2_val ) while XH2O_fluid > X_fluid + 0.01: # too high refinement 1 CO2_val += 0.01 XH2O_fluid = self.get_XH2O_fluid( sample, temperature, pressure, H2O_val, CO2_val ) # ------ Refinement 2 ------ # while XH2O_fluid < X_fluid - 0.001: # too low refinement 2 H2O_val += 0.005 XH2O_fluid = self.get_XH2O_fluid( sample, temperature, pressure, H2O_val, CO2_val ) while XH2O_fluid > X_fluid + 0.001: # too high refinement 2 CO2_val += 0.001 XH2O_fluid = self.get_XH2O_fluid( sample, temperature, pressure, H2O_val, CO2_val ) # ------ Final refinement ------ # while XH2O_fluid < X_fluid - 0.0001: # too low final refinement H2O_val += 0.001 XH2O_fluid = self.get_XH2O_fluid( sample, temperature, pressure, H2O_val, CO2_val ) while XH2O_fluid > X_fluid + 0.0001: # too high final refinement CO2_val += 0.0001 XH2O_fluid = self.get_XH2O_fluid( sample, temperature, pressure, H2O_val, CO2_val ) # ------ Get calculated values ------ # _sample.change_composition( {"H2O": H2O_val, "CO2": CO2_val}, units="wtpt_oxides" ) self.melts.set_bulk_composition( _sample.get_composition(units="wtpt_oxides", normalization="none") ) with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, pressureMPa, initialize=True) (status, temperature, pressureMPa, xmlout) = output[0] fluid_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Fluid") system_mass = self.melts.get_mass_of_phase(xmlout, phase_name="System") liquid_comp = self.melts.get_composition_of_phase( xmlout, phase_name="Liquid", mode="oxide_wt" ) fluid_comp = self.melts.get_composition_of_phase( xmlout, phase_name="Fluid", mode="component" ) if "H2O" in liquid_comp: H2O_liq = liquid_comp["H2O"] else: H2O_liq = 0 if "CO2" in liquid_comp: CO2_liq = liquid_comp["CO2"] else: CO2_liq = 0 if "Water" in fluid_comp: H2O_fl = fluid_comp["Water"] else: H2O_fl = 0.0 if "Carbon Dioxide" in fluid_comp: CO2_fl = fluid_comp["Carbon Dioxide"] else: CO2_fl = 0.0 XH2O_fluid = H2O_fl if verbose: return { "temperature": temperature, "pressure": pressure, "H2O_liq": H2O_liq, "CO2_liq": CO2_liq, "XH2O_fl": H2O_fl, "XCO2_fl": CO2_fl, "FluidProportion_wt": 100 * fluid_mass / system_mass, } if verbose is False: return {"CO2_liq": CO2_liq, "H2O_liq": H2O_liq}
[docs] def calculate_equilibrium_fluid_comp(self, sample, temperature, pressure, verbose=False, **kwargs): """ Returns H2O and CO2 concentrations in wt% in a fluid in equilibrium with the given sample at the given P/T condition. Parameters ---------- sample: Sample class Magma major element composition. temperature: float or int Temperature, in degrees C. presure: float or int Pressure, in bars. verbose: bool OPTIONAL: Default is False. If set to True, returns H2O and CO2 concentration in the fluid in mole fraction, mass of the fluid in grams, and proportion of fluid in the system in wt%. Returns ------- dict A dictionary of fluid composition in wt% with keys 'H2O' and 'CO2' is returned. """ _sample = self.preprocess_sample(sample) bulk_comp_dict = _sample.get_composition(units="wtpt_oxides") if isinstance(temperature, float) or isinstance(temperature, int): pass else: raise core.InputError("temp must be type float or int") if (isinstance(pressure, float) or isinstance(pressure, int) or isinstance(pressure, np.float64)): pass else: raise core.InputError("presure must be type float or int") # Check if only single volatile species is passed. If so, can skip calculations. if bulk_comp_dict["H2O"] == 0: if bulk_comp_dict["CO2"] == 0: if verbose is False: return {"CO2": 0.0, "H2O": 0.0} if verbose: return { "CO2": 0.0, "H2O": 0.0, "FluidMass_grams": 0.0, "FluidProportion_wt": 0.0, } else: if verbose is False: return {"CO2": 1.0, "H2O": 0.0} else: if bulk_comp_dict["CO2"] == 0: if verbose is False: return {"CO2": 0.0, "H2O": 1.0} pressureMPa = pressure / 10.0 self.melts.set_bulk_composition(bulk_comp_dict) with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, pressureMPa, initialize=True) (status, temperature, pressureMPa, xmlout) = output[0] fluid_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Fluid") flsystem_wtper = ( 100 * fluid_mass / (fluid_mass + self.melts.get_mass_of_phase(xmlout, phase_name="Liquid")) ) if fluid_mass > 0.0: fluid_comp = self.melts.get_composition_of_phase( xmlout, phase_name="Fluid", mode="component" ) fluid_comp_H2O = fluid_comp["Water"] fluid_comp_CO2 = fluid_comp["Carbon Dioxide"] else: fluid_comp_H2O = 0 fluid_comp_CO2 = 0 self.melts.set_bulk_composition(self.bulk_comp_orig) # reset if verbose is False: simple_return = {"CO2": fluid_comp_CO2, "H2O": fluid_comp_H2O} simple_return = {key: np.float64(value) for key, value in simple_return.items()} return simple_return if verbose: verbose_return = { "CO2": fluid_comp_CO2, "H2O": fluid_comp_H2O, "FluidMass_grams": fluid_mass, "FluidProportion_wt": flsystem_wtper, } verbose_return = {key: np.float64(value) for key, value in verbose_return.items()} return verbose_return
[docs] def calculate_saturation_pressure( self, sample, temperature, verbose=False, **kwargs ): """ Calculates the saturation pressure of a sample composition. Parameters ---------- sample: Sample class Magma major element composition. temperature: flaot or int Temperature of the sample in degrees C. verbose: bool OPTIONAL: Default is False. If set to False, only the saturation pressure is returned. If set to True, the saturation pressure, mass of fluid in grams, proportion of fluid in wt%, and H2O and CO2 concentrations in the fluid in mole fraction are all returned in a dict. Returns ------- float or dict If verbose is set to False: Saturation pressure in bars. If verbose is set to True: dict of all calculated values. """ _sample = self.preprocess_sample(sample) bulk_comp = _sample.get_composition(units="wtpt_oxides") self.melts.set_bulk_composition(bulk_comp) # Coarse search # NOTE that pressure is in MPa for MagmaSat calculations but reported in bars. pressureMPa = 2000 # Check if saturated at 2000 MPa (rare, for deep samples) with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, pressureMPa, initialize=True) (status, temperature, pressureMPa, xmlout) = output[0] fluid_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Fluid") if fluid_mass <= 0: # if not sat'd at 2000 MPa while fluid_mass <= 0: pressureMPa -= 100 if pressureMPa <= 0: break # composition needs to be reset for each refinement self.melts.set_bulk_composition(bulk_comp) with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, pressureMPa, initialize=True) (status, temperature, pressureMPa, xmlout) = output[0] fluid_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Fluid") fluid_mass = 0 pressureMPa += 100 elif fluid_mass > 0: # if sat'd at 2000 MPa, add pressure while fluid_mass > 0: pressureMPa += 100 # composition needs to be reset for each refinement self.melts.set_bulk_composition(bulk_comp) with redirect_stdout(_f): output = self.melts.equilibrate_tp( temperature, pressureMPa, initialize=True) (status, temperature, pressureMPa, xmlout) = output[0] fluid_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Fluid") fluid_mass = 1.0 pressureMPa -= 100 # Refined search 1 if fluid_mass <= 0: # proceed down pressure search while fluid_mass <= 0: pressureMPa -= 10 if pressureMPa <= 0: break # composition needs to be reset for each refinement self.melts.set_bulk_composition(bulk_comp) with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, pressureMPa, initialize=True) (status, temperature, pressureMPa, xmlout) = output[0] fluid_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Fluid") fluid_mass = 0 pressureMPa += 10 elif fluid_mass > 0: # proceed upward pressure search while fluid_mass > 0: pressureMPa += 10 # composition needs to be reset for each refinement self.melts.set_bulk_composition(bulk_comp) with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, pressureMPa, initialize=True) (status, temperature, pressureMPa, xmlout) = output[0] fluid_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Fluid") fluid_mass = 1.0 pressureMPa -= 10 # Refined search 2 if fluid_mass <= 0: # proceed down pressure search while fluid_mass <= 0: pressureMPa -= 1 if pressureMPa <= 0: break # composition needs to be reset for each refinement self.melts.set_bulk_composition(bulk_comp) with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, pressureMPa, initialize=True) (status, temperature, pressureMPa, xmlout) = output[0] fluid_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Fluid") pass elif fluid_mass > 0: # proceed upward pressure search while fluid_mass > 0: pressureMPa += 1 # composition needs to be reset for each refinement self.melts.set_bulk_composition(bulk_comp) with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, pressureMPa, initialize=True) (status, temperature, pressureMPa, xmlout) = output[0] fluid_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Fluid") if pressureMPa != np.nan: satP = pressureMPa * 10 # convert pressure to bars flmass = fluid_mass flsystem_wtper = (100 * fluid_mass / (fluid_mass + self.melts.get_mass_of_phase(xmlout, phase_name="Liquid"))) flcomp = self.melts.get_composition_of_phase(xmlout, phase_name="Fluid", mode="component") try: flH2O = flcomp["Water"] except Exception: flH2O = 0.0 try: flCO2 = flcomp["Carbon Dioxide"] except Exception: flCO2 = 0.0 else: flmass = np.nan flsystem_wtper = np.nan flH2O = np.nan flCO2 = np.nan warnmessage = "Calculation failed." self.melts.set_bulk_composition( self.bulk_comp_orig ) # this needs to be reset always! if verbose is False: try: w.warn(warnmessage) except Exception: pass return np.float64(satP) elif verbose: try: w.warn(warnmessage) except Exception: pass verbose_return = { "SaturationP_bars": satP, "FluidMass_grams": flmass, "FluidProportion_wt": flsystem_wtper, "XH2O_fl": flH2O, "XCO2_fl": flCO2, } verbose_return = {key: np.float64(value) for key, value in verbose_return.items()} return verbose_return
[docs] def calculate_isobars_and_isopleths( self, sample, temperature, pressure_list, isopleth_list=None, smooth_isobars=True, smooth_isopleths=True, print_status=True, **kwargs ): """ Calculates isobars and isopleths at a constant temperature for a given sample. Isobars can be calculated for any number of pressures. Isobars are calculated using 5 XH2O values (0, 0.25, 0.5, 0.75, 1). Parameters ---------- sample: Sample class Magma major element composition. temperature: float Temperature in degrees C. pressure_list: list or float List of all pressure values at which to calculate isobars, in bars. If only one value is passed it can be as float instead of list. isopleth_list: list or float OPTIONAL: Default value is None in which case only isobars will be calculated. List of all fluid compositions in mole fraction H2O (XH2Ofluid) at which to calcualte isopleths. Values can range from 0-1. If only one value is passed it can be as float instead of list. smooth_isobars: bool OPTIONAL. Default is True. If set to True, polynomials will be fit to the computed isobar points. smooth_isopleths: bool OPTIONAL. Default is True. If set to True, polynomials will be fit to the computed isopleth points. print_status: bool OPTIONAL: Default is True. If set to True, progress of the calculations will be printed to the terminal. Returns ------- pandas DataFrame objects Two pandas DataFrames are returned; the first has isobar data, and the second has isopleth data. Columns in the isobar dataframe are 'Pressure', 'H2Omelt', and 'CO2melt', correpsonding to pressure in bars and dissolved H2O and CO2 in the liquid in wt%. Columns in the isopleth dataframe are 'Pressure', 'H2Ofl', and 'CO2fl', corresponding to pressure in bars and H2O and CO2 concentration in the H2O-CO2 fluid, in wt%. """ if isinstance(pressure_list, list): P_vals = pressure_list elif isinstance(pressure_list, int) or isinstance(pressure_list, float): P_vals = [pressure_list] else: raise core.InputError("pressure_list must be a single float (1000.0), int (1000), or " "list of those [1000, 2000.0, 3000].") if isopleth_list is None: pass elif isinstance(isopleth_list, list): iso_vals = isopleth_list else: iso_vals = [isopleth_list] _sample = self.preprocess_sample(sample) required_iso_vals = [0, 0.25, 0.5, 0.75, 1] all_iso_vals = iso_vals + required_iso_vals all_iso_vals = list(dict.fromkeys(all_iso_vals)) # remove duplicates all_iso_vals.sort() # sort from smallest to largest isobar_data = [] isopleth_data = [] for X in iso_vals: isopleth_data.append([X, 0.0, 0.0]) # Calculate equilibrium phase assemblage for all P/T conditions, check if saturated in # fluid... for i in P_vals: guess = 0.0 if print_status: print("Calculating isobar at " + str(i) + " bars") X_iter = 0 for X in all_iso_vals: X_iter += 1 if print_status: if isopleth_list is not None and X in iso_vals: sys.stdout.write( "\r Calculating isopleth at XH2Ofluid = " + str(X) + " " ) if X not in iso_vals: sys.stdout.write( "\r Calculating isobar control point at XH2Ofluid = " + str(X) + " " ) if X_iter == len(all_iso_vals): sys.stdout.write( "\r done. " " " " \n" ) saturated_vols = self.calculate_dissolved_volatiles( sample=_sample, temperature=temperature, pressure=i, H2O_guess=guess, X_fluid=X, ) if X in required_iso_vals: isobar_data.append([i, saturated_vols["H2O_liq"], saturated_vols["CO2_liq"]]) if X in iso_vals: isopleth_data.append([X, saturated_vols["H2O_liq"], saturated_vols["CO2_liq"]]) guess = saturated_vols["H2O_liq"] if print_status: print("Done!") isobars_df = pd.DataFrame(isobar_data, columns=["Pressure", "H2O_liq", "CO2_liq"]) isopleths_df = pd.DataFrame(isopleth_data, columns=["XH2O_fl", "H2O_liq", "CO2_liq"]) self.melts.set_bulk_composition(self.bulk_comp_orig) # reset if smooth_isobars: isobars_smoothed = vplot.smooth_isobars_and_isopleths(isobars=isobars_df) res_isobars = isobars_smoothed.copy() else: res_isobars = isobars_df.copy() if smooth_isopleths: isopleths_smoothed = vplot.smooth_isobars_and_isopleths(isopleths=isopleths_df) res_isopleths = isopleths_smoothed.copy() else: res_isopleths = isopleths_df.copy() return res_isobars, res_isopleths
[docs] def calculate_degassing_path(self, sample, temperature, pressure="saturation", fractionate_vapor=0.0, init_vapor=0.0, final_pressure=1.0, steps=50, **kwargs): """ Calculates degassing path for one sample Parameters ---------- sample: Sample class Magma major element composition. Legacy info follows, might still be useful? If pulling from an uploaded file with data for many samples, first call get_sample_composition() to get the sample desired. Then pass the result into this function. temperature: float Temperature at which to calculate degassing paths, in degrees C. pressure: string, float, int, list, or numpy array OPTIONAL. The perssure at which to begin the degassing calculations. Default value is 'saturation', which runs the calculation with the initial pressure at the saturation pressure. If a pressure greater than the saturation pressure is input, the calculation will start at saturation, since this is the first pressure at which any degassing will occur. fractionate_vapor: float OPTIONAL. Proportion of vapor removed at each pressure step. Default value is 0.0 (completely closed-system degassing). Specifies the type of calculation performed, either closed system (0.0) or open system (1.0) degassing. If any value between <1.0 is chosen, user can also specify the 'init_vapor' argument (see below). A value in between 0 and 1 will remove that proportion of vapor at each step. For example, for a value of 0.2, the calculation will remove 20% of the vapor and retain 80% of the vapor at each pressure step. init_vapor: float OPTIONAL. Default value is 0.0. Specifies the amount of vapor (in wt%) coexisting with the melt before degassing. final_pressure: float OPTIONAL. The final pressure on the degassing path, in bars. Ignored if a list or numpy array is passed as the pressure variable. Default is 1 bar. steps: int OPTIONAL. Default value is 50. Specifies the number of steps in pressure space at which dissolved volatile concentrations are calculated. Returns ------- pandas DataFrame object """ sys.stdout.write("Finding saturation point... ") # print start of calculation to terminal _sample = self.preprocess_sample(sample) # Normalize sample composition _normed_comp = _sample.get_composition(normalization="standard") _sample.change_composition(_normed_comp) _sample_dict = _sample.get_composition() # ------ RESET MELTS ------ # # MELTS needs to be reloaded here. If an unfeasible composition gets set inside of MELTS, # which can happen when running open-system degassing path calcs, the following calls to # MELTS will fail. This prevents that from happening. self.melts = equilibrate.MELTSmodel("1.2.0") # Suppress phases not required in the melts simulation phases = self.melts.get_phase_names() for phase in phases: self.melts.set_phase_inclusion_status({phase: False}) self.melts.set_phase_inclusion_status({"Fluid": True, "Liquid": True}) self.melts.set_bulk_composition(_sample_dict) # ------------------------- # # Get saturation pressure data = self.calculate_saturation_pressure(sample=_sample, temperature=temperature, verbose=True) if (type(pressure) == str or type(pressure) == float or type(pressure) == int): if pressure == "saturation" or pressure >= data["SaturationP_bars"]: SatP_bars = data["SaturationP_bars"] else: SatP_bars = pressure step_size = (SatP_bars - final_pressure) / steps P_array_bars = np.arange(SatP_bars, final_pressure, -step_size) P_array_MPa = P_array_bars/10.0 # add last few MPa steps P_array_MPa = np.append(P_array_MPa, 0.5) P_array_MPa = np.append(P_array_MPa, 0.1) elif type(pressure) == list: SatP_bars = max(pressure) finalP_bars = min(pressure) step_size = (SatP_bars - finalP_bars) / steps P_array_bars = np.arange(SatP_bars, finalP_bars, -step_size) P_array_MPa = P_array_bars/10.0 elif type(pressure) == np.ndarray: P_array_bars = pressure P_array_MPa = P_array_bars/10.0 # # convert number of steps to step size # MPa_step = SatP_MPa / steps # if MPa_step < 1: # MPa_step = 1 fl_wtper = data["FluidProportion_wt"] while fl_wtper <= init_vapor: with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, np.amax(P_array_MPa), initialize=True) (status, temperature, p, xmlout) = output[0] fl_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Fluid") liq_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Liquid") fl_comp = self.melts.get_composition_of_phase(xmlout, phase_name="Fluid") fl_wtper = 100 * fl_mass / (fl_mass + liq_mass) try: _sample_dict["H2O"] += fl_comp["H2O"] * 0.0005 except Exception: _sample_dict["H2O"] = _sample_dict["H2O"] * 1.1 try: _sample_dict["CO2"] += fl_comp["CO2"] * 0.0005 except Exception: _sample_dict["CO2"] = _sample_dict["CO2"] * 1.1 _sample = sample_class.Sample(_sample_dict) _sample_dict = _sample.get_composition(normalization="standard", units="wtpt_oxides") self.melts.set_bulk_composition(_sample_dict) # reset MELTS press = [] H2Oliq = [] CO2liq = [] H2Ofl = [] CO2fl = [] fluid_wtper = [] iterno = 0 sys.stdout.write("\r") # carriage return to remove previous printed text for i in P_array_MPa: # Handle status_bar iterno += 1 percent = iterno / len(P_array_MPa) batchfile.status_bar.status_bar(percent, btext="Calculating degassing path...") fl_mass = 0.0 self.melts.set_bulk_composition(_sample_dict) with redirect_stdout(_f): output = self.melts.equilibrate_tp(temperature, i, initialize=True) (status, temp, p, xmlout) = output[0] liq_comp = self.melts.get_composition_of_phase(xmlout, phase_name="Liquid") fl_comp = self.melts.get_composition_of_phase(xmlout, phase_name="Fluid", mode="component") liq_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Liquid") fl_mass = self.melts.get_mass_of_phase(xmlout, phase_name="Fluid") fl_wtper = 100 * fl_mass / (fl_mass + liq_mass) if fl_mass > 0: press.append(p * 10.0) try: H2Oliq.append(liq_comp["H2O"]) except Exception: H2Oliq.append(0) try: CO2liq.append(liq_comp["CO2"]) except Exception: CO2liq.append(0) try: H2Ofl.append(fl_comp["Water"]) except Exception: H2Ofl.append(0) try: CO2fl.append(fl_comp["Carbon Dioxide"]) except Exception: CO2fl.append(0) fluid_wtper.append(fl_wtper) try: _sample_dict["H2O"] = (liq_comp["H2O"] + (_sample_dict["H2O"] - liq_comp["H2O"]) * (1.0 - fractionate_vapor)) except Exception: _sample_dict["H2O"] = 0 try: _sample_dict["CO2"] = (liq_comp["CO2"] + (_sample_dict["CO2"] - liq_comp["CO2"]) * (1.0 - fractionate_vapor)) except Exception: _sample_dict["CO2"] = 0 _sample = sample_class.Sample(_sample_dict) _sample_dict = _sample.get_composition(normalization="standard", units="wtpt_oxides") self.melts.set_bulk_composition(self.bulk_comp_orig) # this needs to be reset always! open_degassing_df = pd.DataFrame(list(zip(press, H2Oliq, CO2liq, H2Ofl, CO2fl, fluid_wtper)), columns=["Pressure_bars", "H2O_liq", "CO2_liq", "XH2O_fl", "XCO2_fl", "FluidProportion_wt"]) open_degassing_df = open_degassing_df[open_degassing_df.CO2_liq >= 0.0] open_degassing_df = open_degassing_df[open_degassing_df.H2O_liq >= 0.0] return open_degassing_df