diff --git a/input_files/output.yaml b/input_files/output.yaml index 817360f..24d6e51 100644 --- a/input_files/output.yaml +++ b/input_files/output.yaml @@ -1,8 +1,8 @@ # Graphical output selection -Plt_melt_species: True -Plt_gas_species_wt: False -Plt_gas_species_mol: True +plot_melt_species: True +plot_gas_species_wt: False +plot_gas_species_mol: True -Plt_gas_fraction: True -Plt_fo2_dFMQ: False +plot_gas_fraction: True +plot_fo2_dFMQ: False diff --git a/src/evo/__main__.py b/src/evo/__main__.py index 0b9df21..f9a0939 100644 --- a/src/evo/__main__.py +++ b/src/evo/__main__.py @@ -21,7 +21,10 @@ def main(argv=None): ) my_parser.add_argument( - "-o", "--output", help="the folder location to write the results to" + "-o", + "--output", + help="the folder location to write the results to", + default="outputs", ) # Parse in files diff --git a/src/evo/conversions.py b/src/evo/conversions.py index fd07ecd..6700400 100644 --- a/src/evo/conversions.py +++ b/src/evo/conversions.py @@ -16,10 +16,11 @@ def K2C(T): return T - 273.15 -def truncate(n, decimals=0): - """Round a number mathematically""" - multiplier = 10**decimals - return np.floor(n * multiplier) / multiplier +def round_down(n, dp): + """Round a number according to the number of decimal places""" + step_size = min(1, dp) + factor = 1 / step_size + return np.floor(n * factor) / factor # ------------------------------------------------------------------------ diff --git a/src/evo/dgs.py b/src/evo/dgs.py index 99ca6d6..76494d7 100644 --- a/src/evo/dgs.py +++ b/src/evo/dgs.py @@ -208,6 +208,6 @@ def run_evo(f_chem, f_env, f_out, folder="outputs"): print("Run time is ", end - start) df = writeout_file(sys, gas, melt, sys.P_track) - writeout_figs(sys, melt, gas, out, sys.P_track) + writeout_figs(df, run.results_folder, out=out) return df diff --git a/src/evo/fixed_weights.py b/src/evo/fixed_weights.py index 8fb075c..001fd72 100644 --- a/src/evo/fixed_weights.py +++ b/src/evo/fixed_weights.py @@ -1142,7 +1142,7 @@ def fixed_weights_cohsn(guesses: list[float]): ls.append(np.float64(0)) # Round P down to nearest bar ready for solver to find WgT at 'starting' pressure. - sys.P = cnvs.truncate(sys.P, decimals=0) + sys.P = cnvs.round_down(sys.P, sys.P_step) run.P_START = sys.P return P_sat, values, gamma, mols, melt.graphite_sat diff --git a/src/evo/plots.py b/src/evo/plots.py new file mode 100644 index 0000000..b1e475e --- /dev/null +++ b/src/evo/plots.py @@ -0,0 +1,146 @@ +import matplotlib.pyplot as plt +import numpy as np + + +# fO2 and dFMQ +def plot_fo2FMQ(df): + """Plots the fO2 relative to FMQ against pressure""" + fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True) + ax1.plot(df["FMQ"], df["P"]) + ax1.set_yscale("log") + ax1.set_ylabel("Pressure (bars)") + ax1.invert_yaxis() + ax1.set_xlabel(r"$\Delta$ FMQ") + + fo2 = [] + for x in df["fo2"]: + fo2.append(np.log10(np.exp(x))) + + ax2.plot(fo2, df["P"]) + ax2.set_xlabel("log(10) fO2") + + return fig + + +# Gas speciation mol +def plot_gasspecies_mol(df): + """Plots the gas speciation (mole fraction) vs pressure""" + + fig, ax = plt.subplots() + ax.plot(df["mH2O"], df["P"], c="darkblue", label="H2O") + ax.plot(df["mH2"], df["P"], c="steelblue", label="H2") + ax.plot(df["mO2"], df["P"], c="firebrick", label="O2") + ax.annotate("H2O", xy=[df["mH2O"].iloc[-1], df["P"].iloc[-1]], color="darkblue") + ax.annotate("H2", xy=[df["mH2"].iloc[-1], df["P"].iloc[-1]], color="steelblue") + ax.annotate("O2", xy=[df["mO2"].iloc[-1], df["P"].iloc[-1]], color="firebrick") + + ax.plot(df["mCO2"], df["P"], c="saddlebrown", label="CO2") + ax.plot(df["mCO2"], df["P"], c="darkgreen", label="CO") + ax.plot(df["mCH4"], df["P"], c="forestgreen", label="CH4") + ax.annotate("CO2", xy=[df["mCO2"].iloc[-1], df["P"].iloc[-1]], color="saddlebrown") + ax.annotate("CO", xy=[df["mCO2"].iloc[-1], df["P"].iloc[-1]], color="darkgreen") + ax.annotate("CH4", xy=[df["mCH4"].iloc[-1], df["P"].iloc[-1]], color="forestgreen") + + ax.plot(df["mSO2"], df["P"], c="maroon", label="SO2") + ax.plot(df["mS2"], df["P"], c="goldenrod", label="S2") + ax.plot(df["mH2S"], df["P"], c="pink", label="H2S") + ax.annotate("SO2", xy=[df["mSO2"].iloc[-1], df["P"].iloc[-1]], color="maroon") + ax.annotate("S2", xy=[df["mS2"].iloc[-1], df["P"].iloc[-1]], color="goldenrod") + ax.annotate("H2S", xy=[df["mH2S"].iloc[-1], df["P"].iloc[-1]], color="pink") + + ax.plot(df["mN2"], df["P"], c="grey", label="N2") + ax.annotate("N2", xy=[df["mN2"].iloc[-1], df["P"].iloc[-1]], color="grey") + + ax.set_xscale("log") + ax.invert_yaxis() + ax.set_xlabel("Gas speciation (mol frac)") + ax.set_ylabel("Pressure (bar)") + return fig + + +# Gas speciation wt% + + +def plot_gasspecies_wt(df): + """Plots the gas speciation (weight fraction) vs pressure""" + + fig, ax = plt.subplots() + + for g, color in zip( + ["wH2O", "wH2", "wO2", "wCO2", "wCO", "wCH4", "wSO2", "wS2", "wH2S", "wN2"], + [ + "darkblue", + "steelblue", + "firebrick", + "saddlebrown", + "darkgreen", + "mediumseagreen", + "maroon", + "goldenrod", + "pink", + "grey", + ], + ): + ax.plot(df[g], df["P"], c=color, label=g[1:]) + ax.annotate(g[1:], xy=[df[g].iloc[-1], df["P"].iloc[-1]], color=color) + + ax.set_xscale("log") + ax.invert_yaxis() + ax.set_xlabel("Gas phase (wt %)") + ax.set_ylabel("Pressure (bar)") + return fig + + +# Melt volatile wt% + + +def plot_meltspecies(df): + """Plots the volatile content of the melt vs pressure""" + + fig, ax = plt.subplots() + for s, color in zip( + [ + "H2O_melt", + "H2_melt", + "CO2_melt", + "CO_melt", + "CH4_melt", + "S2-_melt", + "S6+_melt", + "N_melt", + ], + [ + "darkblue", + "steelblue", + "saddlebrown", + "darkgreen", + "mediumseagreen", + "maroon", + "pink", + "grey", + ], + ): + ax.plot(df[s], df["P"], c=color, label=s.split("_")[0]) + ax.annotate(s.split("_")[0], xy=[df[s].iloc[-1], df["P"].iloc[-1]], color=color) + + ax.set_xscale("log") + ax.invert_yaxis() + ax.set_xlabel("Melt volatile content (wt%)") + ax.set_ylabel("Pressure (bar)") + return fig + + +# Exsolved gas mass and volume + + +def plot_gasfraction(df): + """Plots the gas volume fraction vs pressure""" + + fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True) + ax1.plot(df["Gas_wt"], df["P"]) + ax1.invert_yaxis() + ax1.set_xlabel("Exsolved gas wt%") + ax1.set_ylabel("Pressure (bar)") + ax2.plot(df["Exsol_vol%"], df["P"]) + ax2.set_xlabel("Exsolved gas volume %") + return fig diff --git a/src/evo/sat_pressure.py b/src/evo/sat_pressure.py index a0dd835..1d30efe 100644 --- a/src/evo/sat_pressure.py +++ b/src/evo/sat_pressure.py @@ -855,7 +855,7 @@ def quadratic(): ls.append(np.float64(0)) # Round P down to nearest bar ready for solver to find WgT at 'starting' pressure. - sys.P = cnvs.truncate(sys.P, decimals=0) + sys.P = cnvs.round_down(sys.P, sys.P_step) run.P_START = sys.P return P_sat, values, gamma, mols, melt.graphite_sat diff --git a/src/evo/writeout.py b/src/evo/writeout.py index 37e8596..354dee0 100644 --- a/src/evo/writeout.py +++ b/src/evo/writeout.py @@ -8,6 +8,7 @@ import pandas as pd import evo.conversions as cnvt +import evo.plots as plots import evo.sat_pressure as sat @@ -202,7 +203,7 @@ def writeout_file(sys, gas, melt, P, crashed=False): return df -def writeout_figs(sys, melt, gas, out, P): +def writeout_figs(df, filepath, out=None): """ Selects the figures to produce at the end of a successful run. @@ -211,212 +212,36 @@ def writeout_figs(sys, melt, gas, out, P): Parameters ---------- - sys : ThermoSystem class - Active instance of the ThermoSystem class - melt : Melt class - Active instance of the Melt class - gas : Gas class - Active instance of the Gas class + df: + DataFrame containing the output data + filepath : Path + Path to the desired output folder out : Output class or None Active instance of the Output class, or None - P : list of float - List of the pressures each step was calculated at (bar) """ - filepath = sys.run.results_folder - # filelist = glob.glob(str(filepath / "*.png")) filelist = list(filepath.glob("*.png")) # Removes previous files so if output specification is changed # there is no confusion as to up to date files. for file in filelist: file.unlink() - if ( - out is not None - ): # If an output file listing requested figures has been included: - if out.Plt_melt_species: - plot_meltspecies(melt, P, filepath) - - if out.Plt_gas_species_wt: - plot_gasspecies_wt(gas, P, filepath) - - if out.Plt_gas_species_mol: - plot_gasspecies_mol(gas, P, filepath) - - if out.Plt_gas_fraction: - plot_gasfraction(sys, P, filepath) - - if out.Plt_fo2_dFMQ: - plot_fo2FMQ(melt, gas, P, filepath) - - if out.Plt_gas_fraction is not None: - pass - # plot - - else: # else plot every option - plot_meltspecies(melt, P, filepath) - plot_gasspecies_wt(gas, P, filepath) - plot_gasspecies_mol(gas, P, filepath) - plot_gasfraction(sys, P, filepath) - plot_fo2FMQ(melt, gas, P, filepath) - - -# fO2 and dFMQ -def plot_fo2FMQ(melt, gas, P, path): - """Plots the fO2 relative to FMQ against pressure""" - plt.plot(P, melt.fmq) - plt.xscale("log") - plt.xlabel("Pressure (bars)") - plt.ylabel(r"$\Delta$ FMQ") - plt.savefig(path / "FMQ.png") - plt.close() - - fo2 = [] - for x in gas.fo2: - fo2.append(np.log10(np.exp(x))) - - plt.plot(P, fo2) - # plt.xscale('log') - plt.xlabel("Pressure (bars)") - plt.ylabel("log(10) fO2") - plt.savefig(path / "fO2.png") - plt.close() - - -# Gas speciation mol -def plot_gasspecies_mol(gas, P, path): - """Plots the gas speciation (mole fraction) vs pressure""" - - plt.plot(gas.mH2O[1:], P, c="darkblue") - plt.plot(gas.mH2[1:], P, c="steelblue") - plt.plot(gas.mO2[1:], P, c="firebrick") - plt.annotate("H2O", xy=[gas.mH2O[-1], P[-1]], color="darkblue") - plt.annotate("H2", xy=[gas.mH2[-1], P[-1]], color="steelblue") - plt.annotate("O2", xy=[gas.mO2[-1], P[-1]], color="firebrick") - if ( - gas.sys.run.GAS_SYS == "COH" - or gas.sys.run.GAS_SYS == "COHS" - or gas.sys.run.GAS_SYS == "COHSN" - ): - plt.plot(gas.mCO2[1:], P, c="saddlebrown") - plt.plot(gas.mCO[1:], P, c="darkgreen") - plt.plot(gas.mCH4[1:], P, c="forestgreen") - plt.annotate("CO2", xy=[gas.mCO2[-1], P[-1]], color="saddlebrown") - plt.annotate("CO", xy=[gas.mCO[-1], P[-1]], color="darkgreen") - plt.annotate("CH4", xy=[gas.mCH4[-1], P[-1]], color="forestgreen") - if ( - gas.sys.run.GAS_SYS == "SOH" - or gas.sys.run.GAS_SYS == "COHS" - or gas.sys.run.GAS_SYS == "COHSN" - ): - plt.plot(gas.mSO2[1:], P, c="maroon") - plt.plot(gas.mS2[1:], P, c="goldenrod") - plt.plot(gas.mH2S[1:], P, c="pink") - plt.annotate("SO2", xy=[gas.mSO2[-1], P[-1]], color="maroon") - plt.annotate("S2", xy=[gas.mS2[-1], P[-1]], color="goldenrod") - plt.annotate("H2S", xy=[gas.mH2S[-1], P[-1]], color="pink") - if gas.sys.run.GAS_SYS == "COHSN": - plt.plot(gas.mN2[1:], P, c="grey") - plt.annotate("N2", xy=[gas.mN2[-1], P[-1]], color="grey") - plt.xscale("log") - plt.gca().invert_yaxis() - plt.xlabel(f"Speciation in a {gas.sys.run.GAS_SYS} gas (mol frac)") - plt.ylabel("Pressure (bar)") - plt.savefig(path / "speciation(mol).png") - plt.close() - - -# Gas speciation wt% - - -def plot_gasspecies_wt(gas, P, path): - """Plots the gas speciation (weight fraction) vs pressure""" - - plt.plot(gas.Wt["H2O"], P, c="darkblue") - plt.plot(gas.Wt["H2"], P, c="steelblue") - plt.plot(gas.Wt["O2"], P, c="firebrick") - plt.annotate("H2O", xy=[gas.Wt["H2O"][-1], P[-1]], color="darkblue") - plt.annotate("H2", xy=[gas.Wt["H2"][-1], P[-1]], color="steelblue") - plt.annotate("O2", xy=[gas.Wt["O2"][-1], P[-1]], color="firebrick") - if ( - gas.sys.run.GAS_SYS == "COH" - or gas.sys.run.GAS_SYS == "COHS" - or gas.sys.run.GAS_SYS == "COHSN" - ): - plt.plot(gas.Wt["CO2"], P, c="saddlebrown") - plt.plot(gas.Wt["CO"], P, c="darkgreen") - # plt.plot(gas.Wt['CH4'], P, c='mediumseagreen') - plt.annotate("CO2", xy=[gas.Wt["CO2"][-1], P[-1]], color="saddlebrown") - plt.annotate("CO", xy=[gas.Wt["CO"][-1], P[-1]], color="darkgreen") - # plt.annotate('CH4', xy=[gas.mCH4[-1], 0.9], color='mediumseagreen') - if ( - gas.sys.run.GAS_SYS == "SOH" - or gas.sys.run.GAS_SYS == "COHS" - or gas.sys.run.GAS_SYS == "COHSN" - ): - plt.plot(gas.Wt["SO2"], P, c="maroon") - plt.plot(gas.Wt["S2"], P, c="goldenrod") - plt.plot(gas.Wt["H2S"], P, c="pink") - plt.annotate("SO2", xy=[gas.Wt["SO2"][-1], P[-1]], color="maroon") - plt.annotate("S2", xy=[gas.Wt["S2"][-1], P[-1]], color="goldenrod") - plt.annotate("H2S", xy=[gas.Wt["H2S"][-1], P[-1]], color="pink") - if gas.sys.run.GAS_SYS == "COHSN": - plt.plot(gas.Wt["N2"], P, c="grey") - plt.annotate("N2", xy=[gas.Wt["N2"][-1], P[-1]], color="grey") - plt.xscale("log") - plt.gca().invert_yaxis() - plt.xlabel(f"Gas phase speciation of a {gas.sys.run.GAS_SYS} system (wt %)") - plt.ylabel("Pressure (bar)") - plt.savefig(path / "speciation(wt).png") - plt.close() - - -# Melt volatile wt% - - -def plot_meltspecies(melt, P, path): - """Plots the volatile content of the melt vs pressure""" - - plt.plot(melt.h2o, P, c="darkblue") - plt.plot(melt.h2, P, c="steelblue") - plt.annotate("H2O", xy=[melt.h2o[-1], P[-1]], color="darkblue") - plt.annotate("H2", xy=[melt.h2[-1], P[-1]], color="steelblue") - if ( - melt.sys.run.GAS_SYS == "COH" - or melt.sys.run.GAS_SYS == "COHS" - or melt.sys.run.GAS_SYS == "COHSN" - ): - plt.plot(melt.co2, P, c="saddlebrown") - plt.annotate("CO2", xy=[melt.co2[-1], P[-1]], color="saddlebrown") - if ( - melt.sys.run.GAS_SYS == "SOH" - or melt.sys.run.GAS_SYS == "COHS" - or melt.sys.run.GAS_SYS == "COHSN" - ): - plt.plot(melt.s, P, c="maroon") - plt.annotate("S", xy=[melt.s[-1], P[-1]], color="maroon") - if melt.sys.run.GAS_SYS == "COHSN": - plt.plot(melt.n, P, c="grey") - plt.annotate("N", xy=[melt.n[-1], P[-1]], color="grey") - plt.xscale("log") - plt.gca().invert_yaxis() - plt.xlabel("Melt volatile content (wt%)") - plt.ylabel("Pressure (bar)") - plt.savefig(path / "meltspecies.png") - plt.close() - - -# Exsolved gas mass and volume + plot_map = { + "plot_melt_species": ("melt_species.png", plots.plot_meltspecies), + "plot_gas_species_wt": ("speciation(wt).png", plots.plot_gasspecies_wt), + "plot_gas_species_mol": ("speciation(mol).png", plots.plot_gasspecies_mol), + "plot_gas_fraction": ("exsolved_gas.png", plots.plot_gasfraction), + "plot_fo2_dFMQ": ("fo2_fmq.png", plots.plot_fo2FMQ), + } -def plot_gasfraction(sys, P, path): - """Plots the gas volume fraction vs pressure""" + chose_plots = { + k: getattr(out, k) if out is not None else True for k in plot_map.keys() + } - fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True) - ax1.plot(cnvt.frac2perc(sys.WgT[1:]), P) - ax1.invert_yaxis() - ax1.set_xlabel("Exsolved gas wt%") - ax1.set_ylabel("Pressure (bar)") - ax2.plot(cnvt.frac2perc(sys.GvF), P) - ax2.set_xlabel("Exsolved gas volume %") - plt.savefig(path / "exsolved_gas.png") + for key, (filename, plot_func) in plot_map.items(): + if not chose_plots[key]: + continue + fig = plot_func(df) # Call the plotting function + fig.savefig(filepath / filename) + plt.close(fig) # Ensure the figure is closed to free memory