From e52123f5e93425e22ffb8c53578088c78d40f50d Mon Sep 17 00:00:00 2001 From: Pip Liggins Date: Mon, 23 Oct 2023 19:23:47 -0700 Subject: [PATCH 1/6] Add O'Neill 2022 sulfate capacity --- env.yaml | 4 +- solubility_laws.py | 111 +++++++++++++++++++++++++++++++++++++++------ 2 files changed, 98 insertions(+), 17 deletions(-) diff --git a/env.yaml b/env.yaml index a2afe85..5ae793e 100644 --- a/env.yaml +++ b/env.yaml @@ -46,9 +46,9 @@ C_MODEL: burguisser2015 CO_MODEL: armstrong2015 # CH4 model, pick either ardia2013, or None (preferred for higher fO2 cases above ~ IW+1) CH4_MODEL: ardia2013 -# S models, pick from: sulfide capacity: oneill2002, oneill2020; sulfate capacity: nash2019 ; SCSS: liu2007 +# S models, pick from: sulfide capacity: oneill2002, oneill2020; sulfate capacity: nash2019 , oneill2022 ; SCSS: liu2007 SULFIDE_CAPACITY: oneill2020 -SULFATE_CAPACITY: nash2019 +SULFATE_CAPACITY: oneill2022 SCSS: liu2007 # N model, pick from: libourel2003 (N2 only) N_MODEL: libourel2003 diff --git a/solubility_laws.py b/solubility_laws.py index e247004..4d4d9f6 100644 --- a/solubility_laws.py +++ b/solubility_laws.py @@ -768,12 +768,12 @@ def oneill2020(T, melt): In: Magma Redox Chemistry. """ if not bool(melt.cm_dry): - comp = melt.Cm() + comp = melt.Cm().copy() else: - comp = melt.cm_dry + comp = melt.cm_dry.copy() if comp["fe2o3"] > 0: - comp["feo"] = comp["feo"] + comp["fe2o3"] * 0.8998 + comp["feo"] = comp["feo"] + comp["fe2o3"] * 2 comp["fe2o3"] = 0.0 comp = cnvs.mol2wt(comp) @@ -814,6 +814,67 @@ def oneill2020(T, melt): return capacity +def oneill2022(T, melt): + """ + Returns the sulfate capacity of a melt. + + Parameters + ---------- + T : float + Temperature (K) + melt : Melt class + Active instance of the Melt class + + Returns + ------- + capacity : float + Sulfate (S6+) capacity of melt as a weight fraction + + References + ---------- + O'Neill, H.S.C., Mavrogenes, J.A (2022) The sulfate capacities of silicate melts. + In: GCA. + """ + + if not bool(melt.cm_dry): + comp = melt.Cm().copy() + else: + comp = melt.cm_dry.copy() + + comp = cnvs.mol2wt(comp) + + # Converts the dry wt fraction to normalized mol fraction on single cation basis. + comp = cnvs.single_cat(comp) + + Na = comp["na2o"] + K = comp["k2o"] + Ca = comp["cao"] + Mg = comp["mgo"] + Al = comp["al2o3"] + Fe2 = comp["feo"] + Mn = comp["mno"] + + capacity = ( + np.exp( + -8.02 + + ( + 21100 + + 44000 * Na + + 18700 * Mg + + 4300 * Al + + 44200 * K + + 35600 * Ca + + 12600 * Mn + + 16500 * Fe2 + ) + / T + ) + / 1000000 + ) # Convert ppm -> wt fraction + + return capacity + + # ------------------------------------------------------------------------------ # MELT CONTENT MODEL SELECTION # ------------------------------------------------------------------------------ @@ -1023,7 +1084,7 @@ def sulfate_melt(fS2, fO2, P, T, melt, run, name="nash2019"): Active instance of the Melt class run : RunDef class Active instance of the RunDef class - name : {'nash2019'} + name : {'nash2019', 'oneill2022'} The name of the solubility law to be used Returns @@ -1037,6 +1098,20 @@ def sulfate_melt(fS2, fO2, P, T, melt, run, name="nash2019"): return ratio * (sulfide_melt(fS2, fO2, P, T, melt, name=run.SULFIDE_CAPACITY)) + elif name == "oneill2022": + C6 = oneill2022(T, melt) + + if run.SULFIDE_CAPACITY == "oneill2002": + C2 = oneill2002(fO2, P, melt) + elif run.SULFIDE_CAPACITY == "oneill2020": + C2 = oneill2020(T, melt) + + K = np.exp(55921 / T - 25.07 + 0.6465 * np.log(T)) + + ratio = (C6 / C2) * K * fO2**2 + + return ratio * (sulfide_melt(fS2, fO2, P, T, melt, name=run.SULFIDE_CAPACITY)) + def sulfide_melt(fS2, fO2, P, T, melt, name="oneill2020"): """ @@ -1287,20 +1362,26 @@ def S2_fugacity( The fugacity of S2 (bar) """ - if sulfatename == "nash2019": - # sets the Fe2/Fe3 ratio and dry melt chemistry for the sulfide capacity - melt.cm_dry = melt.iron_fraction(np.log(fO2), ppa=P * 1e5)[0] + # sets the Fe2/Fe3 ratio and dry melt chemistry for the sulfide capacity + melt.cm_dry = melt.iron_fraction(np.log(fO2), ppa=P * 1e5)[0] - # need to find the ratio of s6+/s2-, then use the amount of s2- to get mS2. + if sulfidename == "oneill2002": + C2 = oneill2002(fO2, P, melt) + + elif sulfidename == "oneill2020": + C2 = oneill2020(T, melt) + + # need to find the ratio of s6+/s2-, then use the amount of s2- to get mS2. + if sulfatename == "nash2019": ratio = nash2019(fO2, P, T, melt, run) # s6+/s2- - s2_melt = s_melt * (1 / (1 + ratio)) - if sulfidename == "oneill2002": - capacity = oneill2002(fO2, P, melt) + elif sulfatename == "oneill2022": + C6 = oneill2022(T, melt) + K = np.exp(55921 / T - 25.07 + 0.6465 * np.log(T)) + ratio = (C6 / C2) * K * fO2**2 - elif sulfidename == "oneill2020": - capacity = oneill2020(T, melt) + s2_melt = s_melt * (1 / (1 + ratio)) - fS2 = (s2_melt / capacity) ** (1 / 0.5) * fO2 + fS2 = (s2_melt / C2) ** (1 / 0.5) * fO2 - return fS2 + return fS2 From 7e592dd9d52b2f8e4a4b0fb80b2d4ab0234692bf Mon Sep 17 00:00:00 2001 From: Pip Liggins Date: Thu, 1 Feb 2024 16:25:35 +0000 Subject: [PATCH 2/6] Allow Fe3/FeT ratio as input. Version used for Fuego October paper outputs. --- conversions.py | 35 ++++++++++++++++++++--- dgs_classes.py | 74 ++++++++++++++++++++++++++++++++++++++++++------- env.yaml | 3 ++ ferric.py | 18 +++++++++--- readin.py | 4 ++- sat_pressure.py | 16 +++++++++-- writeout.py | 1 + 7 files changed, 129 insertions(+), 22 deletions(-) diff --git a/conversions.py b/conversions.py index d10fbc0..45f454e 100644 --- a/conversions.py +++ b/conversions.py @@ -117,7 +117,7 @@ def wt2mol(c, *args): def fo2_2F(cm, t, p, lnfo2, model="kc1991"): """ - Generate the fe3/fe2 ratio based on the set fO2 and fO2 model. + Generate the Fe2O3/FeO ratio based on the set fO2 and fO2 model. Parameters ---------- @@ -146,7 +146,7 @@ def fo2_2F(cm, t, p, lnfo2, model="kc1991"): def c2fo2(cm, t, p, model): """ - Calculates fO2 based on Fe3/Fe2 ratio from the melt major element + Calculates fO2 based on Fe2O3/FeO mole ratio from the melt major element composition and the fO2 model. Parameters @@ -173,13 +173,13 @@ def c2fo2(cm, t, p, model): def single_cat(c): - """Convert dry weight fraction to single cation mol fraction. All iron as FeOt""" + """Convert dry weight fraction to single cation mol fraction.""" mol = {} sm = 0 for x in c: - if x in ["al2o3", "na2o", "k2o"]: + if x in ["al2o3", "na2o", "k2o", "fe2o3"]: mol[x] = c[x] / (cnst.m[x] * 0.5) sm += c[x] / (cnst.m[x] * 0.5) elif x in ["mgo", "sio2", "cao", "tio2", "mno", "feo"]: @@ -336,6 +336,33 @@ def generate_fo2_buffer(sys, fo2, P, buffer_choice=None): return fo2_2nno(np.log10(fo2), sys.T, P, sys.run.FMQ_MODEL) +def generate_fo2_fe3fet(sys, melt, P, Fe3_fet): + """ + Returns fO2 according to the fe3/fet ratio + + Parameters + ---------- + sys : ThermoSystem class + Active instance of the ThermoSystem class + melt : Melt class + Active instance of the Melt class + P : float + Pressure (pascal) + F : float + The mFe2O3/mFeO ratio + + Returns + ------- + float + Absolute fO2 + """ + fe3_fe2 = Fe3_fet / (1 - Fe3_fet) + # F = Xfe2o3/Xfeo + F = fe3_fe2 / 2 + + return ferric.kc91_fe3(melt.Cm(), sys.T, P, F=F) + + # ------------------------------------------------------------------------ # CONVERSIONS FOR THE SOLVER SCRIPT # ------------------------------------------------------------------------ diff --git a/dgs_classes.py b/dgs_classes.py index a9cc5d1..c9e59d2 100644 --- a/dgs_classes.py +++ b/dgs_classes.py @@ -104,6 +104,8 @@ class RunDef: Mineral buffer to cite fO2 relative to FO2_buffer_START : float System fO2 relative to `FO2_buffer`. + FE3_FET: float + Set the redox state with the Fe3/FeT ratio FO2_START : float Set the absolute fO2. Do not use in conjunction with `FO2_buffer_START`. FH2_START : float @@ -173,12 +175,16 @@ def __init__(self): self.SCSS = str("liu2007") self.N_MODEL = str("libourel2003") - # initial gas fugacities + # initial redox self.FO2_buffer_SET = False self.FO2_buffer = "FMQ" self.FO2_buffer_START = 0.0 # delta buffer value + self.FE3_FET_SET = False + self.FE3_FET_START = 0.0 # ratio self.FO2_SET = False self.FO2_START = 0.0 # bar + + # initial gas fugacities self.FH2_SET = False self.FH2_START = 0.24 # bar self.FH2O_SET = False @@ -454,14 +460,16 @@ def norm_with_gas(self, melt): # for iron equilibration. # returns the 9 major oxides as weight %, with feoT and volatiles held constant. oxides = cnvs.norm(C_s, (1 - cons)) - self.Cw.update( - oxides - ) # adds the normalised oxide values to the total system mass frac dictionary + # add the normalised oxide values to the total system mass frac dictionary + self.Cw.update(oxides) melt.cm_dry, F = melt.iron_fraction( self.FO2 ) # this is now a duplicate of the one inside init cons, hopefully. melt.F.append(F) + melt.Fe3FeT.append( + melt.cm_dry["fe2o3"] * 2 / (melt.cm_dry["fe2o3"] * 2 + melt.cm_dry["feo"]) + ) def get_wto_tot(self, melt): """ @@ -561,8 +569,8 @@ def fe_save(self, melt, mo2, O2): fo2 = np.log(O2.Y * mo2 * self.P) - # returns anhydrous silicate melt comp mol fraction, Fe2/Fe3. - F = melt.iron_fraction(fo2)[1] + # returns anhydrous silicate melt comp mol fraction, Fe2O3/FeO. + cm_dry, F = melt.iron_fraction(fo2) ofe = ( cnst.m["o"] @@ -574,6 +582,7 @@ def fe_save(self, melt, mo2, O2): self.atomicM["o"] = self.atomicM["o_tot"] - ofe melt.ofe.append(ofe) melt.F.append(F) + melt.Fe3FeT.append(cm_dry["fe2o3"] * 2 / (cm_dry["fe2o3"] * 2 + cm_dry["feo"])) def mass_conservation_reset(self, melt, gas): """ @@ -620,7 +629,9 @@ def mass_conservation_reset(self, melt, gas): melt.s, melt.n, melt.F, + melt.Fe3FeT, melt.ofe, + melt.rho_store, ] gas_wts_f = ["H2O", "O2", "H2", "CO", "CO2", "CH4", "S2", "SO2", "H2S", "N2"] @@ -937,6 +948,7 @@ def __init__(self, run, sys): self.s = [] self.n = [] # stores the melt n content, independent of speciation. self.F = [] # stores the mfe2o3/mFeO ratio at each step. + self.Fe3FeT = [] # stores the fe3+/fe(total) ratio self.ofe = [] self.graph_current = ( 0 # number of moles of graphite in the system at current timestep @@ -977,9 +989,35 @@ def chem_init(self, eles, chem): sm += tmp_Cw[e] length += 1 - # normalise to 100 and define initial system element mass + # if ( + # self.run.FE3_FET_SET is True + # ): # PL: Rather than do this, maybe store the original values to allow the ratio to be calculated later? It's not quite accurate when calculating the saturation P. + # # set correct fe2o3 and feo + # eles.append("fe2o3") + + # # fe3_fe2 = 1 / (1 / self.run.FE3_FET_START - 1) + # fe3_fe2 = self.run.FE3_FET_START / (1 - self.run.FE3_FET_START) + # F = fe3_fe2 / 2 + + # tmp_Cm = cnvs.wt2mol(tmp_Cw) + # tmp_Cm["fe2o3"] = F * tmp_Cm["feo"] / (2 * F + 1) + # tmp_Cm["feo"] = tmp_Cm["feo"] / (1 + 2 * F) + + # # X * molecular weight + # x_mw = {k: v * cnst.m[k] for k, v in tmp_Cm.items()} + + # # wt% normalised + # for e in eles: + # tmp_Cw[e] = x_mw[e] * 100.0 / sum(x_mw.values()) + # else: + # # normalise to 100 + # for e in eles: + # tmp_Cw[e] = tmp_Cw[e] * 100.0 / sm for e in eles: tmp_Cw[e] = tmp_Cw[e] * 100.0 / sm + + # define initial system element mass + for e in eles: self.C_s[e] = tmp_Cw[e] / 100.0 * self.run.MASS # actual oxide mass # add Fe+++ in if component has not been specified @@ -1012,7 +1050,13 @@ def chem_init(self, eles, chem): # if Fe3 and Fe2 exist, then need to determine system fO2; # also updates the Fe mass in system. - if "feo" in eles and "fe2o3" in eles and self.sys.FO2 is None: + if self.sys.FO2 is None and self.run.FE3_FET_SET: + self.sys.FO2 = np.log( + cnvs.generate_fo2_fe3fet( + self.sys, self, self.sys.Ppa, self.run.FE3_FET_START + ) + ) + elif "feo" in eles and "fe2o3" in eles and self.sys.FO2 is None: self.sys.FO2 = np.log( cnvs.c2fo2(self.Cm(), self.sys.T, self.sys.Ppa, self.run.FO2_MODEL) ) @@ -1052,7 +1096,7 @@ def iron_fraction(self, fo2, ppa=None): mFeoT = composition["feo"] - mfe2o3 = F * mFeoT / (1 + 0.8998 * F) + mfe2o3 = F * mFeoT / (1 + 2 * F) mfeo = mfe2o3 / F @@ -1060,7 +1104,7 @@ def iron_fraction(self, fo2, ppa=None): composition["fe2o3"] = mfe2o3 - composition = cnvs.norm(composition) + composition = cnvs.norm(composition, cons=["feo", "fe2o3"]) return composition, F def Cw(self): @@ -1397,6 +1441,11 @@ def melt_composition(self, gas, mols): if self.run.FE_SYSTEM is False: # otherwise this is done in fe_save self.F.append(F) + self.Fe3FeT.append( + self.cm_dry["fe2o3"] + * 2 + / (self.cm_dry["fe2o3"] * 2 + self.cm_dry["feo"]) + ) self.sulfide.append( sl.sulfide_melt( @@ -1429,6 +1478,11 @@ def melt_composition(self, gas, mols): fo2 = O2.Y * gas.mO2[-1] * self.sys.P self.cm_dry, F = self.iron_fraction(np.log(fo2)) self.F.append(F) + self.Fe3FeT.append( + self.cm_dry["fe2o3"] + * 2 + / (self.cm_dry["fe2o3"] * 2 + self.cm_dry["feo"]) + ) self.sulfide.append(0.0) self.sulfate.append(0.0) diff --git a/env.yaml b/env.yaml index 5ae793e..d8cdbe8 100644 --- a/env.yaml +++ b/env.yaml @@ -59,6 +59,9 @@ FO2_buffer_SET: False FO2_buffer: FMQ FO2_buffer_START: 0 +FE3_FET_SET: False +FE3_FET_START: 0.1 + FO2_SET: False FO2_START: 8.3777516486972E-12 diff --git a/ferric.py b/ferric.py index 203fde6..8d4114d 100644 --- a/ferric.py +++ b/ferric.py @@ -77,7 +77,7 @@ def kc91_fo2(C, T, P, lnfo2): return F -def kc91_fe3(Cm, T, P): +def kc91_fe3(Cm, T, P, F=None): """ Calculates the oxygen fugacity (fO2) of a melt given where the ferric/ ferrous ratio is known. @@ -86,11 +86,13 @@ def kc91_fe3(Cm, T, P): ---------- C : dictionary Major element composition of the silicate melt as mole fractions - Required species: Al2O3, FeO, Fe2O3, CaO, Na2O, K2O + Required species: Al2O3, FeO, (Fe2O3), CaO, Na2O, K2O T : float Temperature in degrees K P : float Pressure in pascals (Pa) + F : float, optional + If provided, the fe2o3/feo ratio to use instead of the composition Returns ------- @@ -121,13 +123,21 @@ def kc91_fe3(Cm, T, P): T0 = 1673.0 # K + if not F: + assert Cm["fe2o3"] != 0.0 + F = Cm["fe2o3"] / Cm["feo"] + Cm_feo = Cm["feo"] + Cm["fe2o3"] * 0.8998 + else: + assert ("fe2o3" not in Cm) or (Cm["fe2o3"] == 0.0) + Cm_feo = Cm["feo"] + FO2 = np.exp( ( - np.log(Cm["fe2o3"] / Cm["feo"]) + np.log(F) - b / T - c - dal2o3 * Cm["al2o3"] - - dfeo * (Cm["feo"] + Cm["fe2o3"] * 0.8998) + - dfeo * (Cm_feo) - dcao * Cm["cao"] - dna2o * Cm["na2o"] - dk2o * Cm["k2o"] diff --git a/readin.py b/readin.py index b7e0e04..b4ca4bf 100644 --- a/readin.py +++ b/readin.py @@ -148,7 +148,9 @@ def readin_chem(f, run, sys): # test for iron/ fO2 definition if "feo" in ele_names and "fe2o3" not in ele_names and run.FO2_SET is not True: - if run.GAS_SYS == "OH" and (run.FH2_SET): + if run.FE3_FET_SET is True: + pass + elif run.GAS_SYS == "OH" and (run.FH2_SET): pass elif run.FH2_SET and (run.FH2O_SET or run.WTH2O_SET): pass diff --git a/sat_pressure.py b/sat_pressure.py index 7b1209f..ef9ddb8 100644 --- a/sat_pressure.py +++ b/sat_pressure.py @@ -239,6 +239,8 @@ def get_f(P, sys, melt, gamma): fo2 = np.exp( cnvs.generate_fo2(sys, sys.run.FO2_buffer_START, sys.run.FO2_buffer, P) ) + elif run.FE3_FET_SET: + fo2 = cnvs.generate_fo2_fe3fet(sys, melt, P * 1e5, run.FE3_FET_START) else: fo2 = cnvs.c2fo2(melt.Cm(), sys.T, P * 1e5, sys.run.FO2_MODEL) @@ -351,6 +353,11 @@ def find_p(P, sys, melt): sys, sys.run.FO2_buffer_START, sys.run.FO2_buffer, sys.P ) + elif run.FE3_FET_SET: + sys.FO2 = np.log( + cnvs.generate_fo2_fe3fet(sys, melt, sys.Ppa, run.FE3_FET_START) + ) + else: sys.FO2 = np.log(cnvs.c2fo2(melt.Cm(), sys.T, sys.Ppa, sys.run.FO2_MODEL)) @@ -1065,9 +1072,9 @@ def satp_writeout(sys, melt, gas, P, values, gamma, mols, graph_sat=False): or sys.run.GAS_SYS == "COHSN" ): fo2 = o2y * mO2 * P - melt.cm_dry, F = melt.iron_fraction( - np.log(fo2), ppa=P * 1e5 - ) # Set the Fe2/Fe3 ratio and dry melt chemistry for the sulfide capacity + # Set the Fe2/Fe3 ratio and dry melt chemistry for the sulfide capacity + cm_dry, F = melt.iron_fraction(np.log(fo2), ppa=P * 1e5) + melt.cm_dry.update(cm_dry) # use mS2 to work back and find amount of S2- in melt s2_melt = ( @@ -1132,6 +1139,9 @@ def satp_writeout(sys, melt, gas, P, values, gamma, mols, graph_sat=False): fo2_buffer: cnvs.generate_fo2_buffer(sys, (o2y * mO2 * P), P), "fo2": fo2, "F": F, + "Fe3FeT": melt.cm_dry["fe2o3"] + * 2 + / (melt.cm_dry["fe2o3"] * 2 + melt.cm_dry["feo"]), "rho_bulk": rho_bulk, "rho_melt": rho_melt, "Exsol_vol%": GvF * 100, diff --git a/writeout.py b/writeout.py index 6138b12..4f090a0 100644 --- a/writeout.py +++ b/writeout.py @@ -37,6 +37,7 @@ def get_data(sys, gas, melt, P): ], "fo2": [np.exp(gas.fo2[i]) for i in range(len(P))], "F": melt.F[1:] if not sys.run.SINGLE_STEP else [melt.F[-1]], + "Fe3FeT": melt.Fe3FeT[1:] if not sys.run.SINGLE_STEP else [melt.Fe3FeT[-1]], "rho_bulk": sys.rho, "rho_melt": melt.rho_store, "Exsol_vol%": [(sys.GvF[i] * 100) for i in range(len(P))], From 8c13d0d556480dc293dec8655806a48b745d82c6 Mon Sep 17 00:00:00 2001 From: Pip Liggins Date: Mon, 13 May 2024 10:43:28 +0100 Subject: [PATCH 3/6] Fix writeout errors for non-COHS systems --- sat_pressure.py | 8 ++++---- writeout.py | 49 +++++++++++++++++++++++++++++++++++++------------ 2 files changed, 41 insertions(+), 16 deletions(-) diff --git a/sat_pressure.py b/sat_pressure.py index ef9ddb8..4fbdc34 100644 --- a/sat_pressure.py +++ b/sat_pressure.py @@ -1176,10 +1176,10 @@ def satp_writeout(sys, melt, gas, P, values, gamma, mols, graph_sat=False): "fH2S": (h2sy * mH2S * P), "fS2": (s2y * mS2 * P), "fN2": (n2y * mN2 * P), - "mCO2/CO": gas.mCO2[0] / gas.mCO[0], - "mCO2/H2O": gas.mCO2[0] / gas.mH2O[0], - "mCO2/SO2": gas.mCO2[0] / gas.mSO2[0], - "mH2S/SO2": gas.mH2S[0] / gas.mSO2[0], + "mCO2/CO": gas.mCO2[0] / gas.mCO[0] if "C" in sys.run.GAS_SYS else np.nan, + "mCO2/H2O": gas.mCO2[0] / gas.mH2O[0] if "C" in sys.run.GAS_SYS else np.nan, + "mCO2/SO2": gas.mCO2[0] / gas.mSO2[0] if "S" in sys.run.GAS_SYS else np.nan, + "mH2S/SO2": gas.mH2S[0] / gas.mSO2[0] if "S" in sys.run.GAS_SYS else np.nan, "H2O_melt": melt_h2o * 100, "H2_melt": melt_h2 * 100, "CO2_melt": melt_co2 * 100, diff --git a/writeout.py b/writeout.py index 4f090a0..248d4b3 100644 --- a/writeout.py +++ b/writeout.py @@ -2,6 +2,7 @@ Writes the results of the run to a file and contains options to produce a graph of the results. """ + import conversions as cnvt import numpy as np import pandas as pd @@ -72,18 +73,42 @@ def get_data(sys, gas, melt, P): "fH2S": gas.f["H2S"], "fS2": gas.f["S2"], "fN2": gas.f["N2"], - "mCO2/CO": [gas.mCO2[i + 1] / gas.mCO[i + 1] for i in range(len(P))] - if not sys.run.SINGLE_STEP - else [gas.mCO2[-1] / gas.mCO[-1]], - "mCO2/H2O": [gas.mCO2[i + 1] / gas.mH2O[i + 1] for i in range(len(P))] - if not sys.run.SINGLE_STEP - else [gas.mCO2[-1] / gas.mH2O[-1]], - "mCO2/SO2": [gas.mCO2[i + 1] / gas.mSO2[i + 1] for i in range(len(P))] - if not sys.run.SINGLE_STEP - else [gas.mCO2[-1] / gas.mSO2[-1]], - "mH2S/SO2": [gas.mH2S[i + 1] / gas.mSO2[i + 1] for i in range(len(P))] - if not sys.run.SINGLE_STEP - else [gas.mH2S[-1] / gas.mSO2[-1]], + "mCO2/CO": ( + np.nan + if "C" not in sys.run.GAS_SYS + else ( + [gas.mCO2[i + 1] / gas.mCO[i + 1] for i in range(len(P))] + if not sys.run.SINGLE_STEP + else [gas.mCO2[-1] / gas.mCO[-1]] + ) + ), + "mCO2/H2O": ( + np.nan + if "C" not in sys.run.GAS_SYS + else ( + [gas.mCO2[i + 1] / gas.mH2O[i + 1] for i in range(len(P))] + if not sys.run.SINGLE_STEP + else [gas.mCO2[-1] / gas.mH2O[-1]] + ) + ), + "mCO2/SO2": ( + np.nan + if "S" not in sys.run.GAS_SYS + else ( + [gas.mCO2[i + 1] / gas.mSO2[i + 1] for i in range(len(P))] + if not sys.run.SINGLE_STEP + else [gas.mCO2[-1] / gas.mSO2[-1]] + ) + ), + "mH2S/SO2": ( + np.nan + if "S" not in sys.run.GAS_SYS + else ( + [gas.mH2S[i + 1] / gas.mSO2[i + 1] for i in range(len(P))] + if not sys.run.SINGLE_STEP + else [gas.mH2S[-1] / gas.mSO2[-1]] + ) + ), "H2O_melt": [melt.h2o[i] * 100 for i in range(len(P))], "H2_melt": [melt.h2[i] * 100 for i in range(len(P))], "CO2_melt": [melt.co2[i] * 100 for i in range(len(P))], From 9aef471e24bc028fa99a240618ff5d0335b87417 Mon Sep 17 00:00:00 2001 From: Pip Liggins Date: Mon, 27 May 2024 14:25:50 +0100 Subject: [PATCH 4/6] Update gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index b2d2e4a..19d39f7 100644 --- a/.gitignore +++ b/.gitignore @@ -18,3 +18,5 @@ __pycache__/ .vscode Output/* + +.DS_Store From 76b23d978c9b2eb502bf7b065cb93297704ad015 Mon Sep 17 00:00:00 2001 From: Pip Liggins Date: Wed, 26 Feb 2025 16:12:23 +0000 Subject: [PATCH 5/6] fix failing tests - don't hold iron constant when normalising molar melt composition --- src/evo/dgs_classes.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/evo/dgs_classes.py b/src/evo/dgs_classes.py index 67becaf..ca0dfec 100644 --- a/src/evo/dgs_classes.py +++ b/src/evo/dgs_classes.py @@ -449,9 +449,8 @@ def norm_with_gas(self, melt): # for iron equilibration. # add the normalised oxide values to the total system mass frac dictionary self.Cw.update(oxides) - melt.cm_dry, F = melt.iron_fraction( - self.FO2 - ) # this is now a duplicate of the one inside init cons, hopefully. + # this is now a duplicate of the one inside init cons, hopefully. + melt.cm_dry, F = melt.iron_fraction(self.FO2) melt.F.append(F) melt.Fe3FeT.append( melt.cm_dry["fe2o3"] * 2 / (melt.cm_dry["fe2o3"] * 2 + melt.cm_dry["feo"]) @@ -1094,7 +1093,9 @@ def iron_fraction(self, fo2, ppa=None): composition["fe2o3"] = mfe2o3 - composition = cnvs.norm(composition, cons=["feo", "fe2o3"]) + # PL: need to think through why you would/wouldn't normalise while holding + # Fe constant... + composition = cnvs.norm(composition) # , cons=["feo", "fe2o3"] return composition, F def Cw(self): From 7e43068918994d2f521b8ef40f2b32fcb94657af Mon Sep 17 00:00:00 2001 From: Pip Liggins Date: Wed, 2 Apr 2025 16:46:00 -0700 Subject: [PATCH 6/6] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index f737503..726662e 100644 --- a/README.md +++ b/README.md @@ -32,7 +32,7 @@ and run using Or EVo can be run directly from the terminal: ``` -python evo input_files/chem.yaml input_files/env.yaml --output-options input_files/output.yaml +python -m evo input_files/chem.yaml input_files/env.yaml --output-options input_files/output.yaml ``` The model should run and produce an output file `outputs/dgs_output_*.csv`, and a set of graphs in an 'outputs' folder, if a decompression run has been selected.