diff --git a/scripts/analyser/analyse_unsaturated_hdf5.py b/scripts/analyser/analyse_unsaturated_hdf5.py index 861c1b9..dc8101c 100644 --- a/scripts/analyser/analyse_unsaturated_hdf5.py +++ b/scripts/analyser/analyse_unsaturated_hdf5.py @@ -1,287 +1,726 @@ # ----------------------------------------------------------------------------- #Copyright (c) 2015 Fabien Georget , Princeton #University #All rights reserved. # #Redistribution and use in source and binary forms, with or without #modification, are permitted provided that the following conditions are met: # #1. Redistributions of source code must retain the above copyright notice, this #list of conditions and the following disclaimer. # #2. Redistributions in binary form must reproduce the above copyright notice, #this list of conditions and the following disclaimer in the documentation #and/or other materials provided with the distribution. # #3. Neither the name of the copyright holder nor the names of its contributors #may be used to endorse or promote products derived from this software without #specific prior written permission. # #THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND #ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED #WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE #DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE #FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL #DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR #SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER #CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, #OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE #OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. # ----------------------------------------------------------------------------- print(" * load modules") import sys import h5py import numpy as np + import matplotlib.pyplot as plt +import scipy.optimize as sco + +#from matplotlib.figure import Figure + + if len(sys.argv) < 2: print("One argument required, aborting") exit(1) def fail(message, exception, err_code=4): print("{0} -> \n{1}".format(message, str(exception))) import os os._exit(err_code) + +def diffusion_profile(sqr_t, a): + return a*sqr_t + + +class UnsaturatedHDF5File: + """Unsaturated driver result file""" + def __init__(self, filename): + self.load_file(filename) + self.parse_timesteps() + self.parse_database() + self.parse_mesh() + self.parse_gas_nodes() + self.parse_ph() + + + def load_file(self, filename): + try: + self.f = h5py.File(filename, mode='r') + self.filename = filename + except Exception as e: + fail("Error while opening file : {0}".format(filename), e) + + + def reload(self): + self.f.close() + try: + self.f = h5py.File(self.filename, mode='r') + except Exception as e: + fail("Error while opening file : {0}".format(self.filename), e) + self.parse_timesteps() + + + def parse_timesteps(self): + try: + self.str_ts = [x for x in self.f][:-2] + self.str_ts.sort(key=float) + self.ts = [float(t) for t in self.str_ts] + self.map_ts = dict(zip(self.ts, self.str_ts)) + except Exception as e: + fail("Error while parsing timesteps", e) + + def get_list_species(self, species): + return list(self.database[species][:]) + + + def parse_database(self): + try: + self.database = self.f['database'] + except Exception as e: + fail("Error while parsing database section", e) + + try: + self.list_components = self.get_list_species('components') + self.list_minerals = self.get_list_species('minerals') + self.list_gas = self.get_list_species('gas') + self.list_aqueous = self.get_list_species('aqueous') + except Exception as e: + fail("Error while parsing database labels", e) + + + def parse_mesh(self): + try: + self.mesh = self.f['mesh'] + self.coordinates = [x for x in self.mesh['coordinates'][:]] + self.nb_nodes = len(self.coordinates) + self.nodes = [int(x) for x in range(self.nb_nodes)] + except Exception as e: + fail("Error while parsing mesh", e) + + + def get_ph_main(self, timestep, node): + if (not self.is_gas_node[node]): + chem_sol = self.get_timestep(timestep)['chemistry_solution'][str(node)] + log_h = chem_sol['main_variables'][self.id_h] + log_gamma_h = chem_sol['log_gamma'][self.id_h] + return -(log_gamma_h + log_h) + else: + return 0; + + + def get_ph_second(self, timestep, node): + if (not self.is_gas_node[node]): + chem_sol = self.get_timestep(timestep)['chemistry_solution'][str(node)] + log_h = np.log10(chem_sol['secondary_molalities'][self.id_h]) + log_gamma_h = chem_sol['log_gamma'][self.id_logh] + return -(log_gamma_h + log_h) + else: + return 0; + + def parse_ph(self): + if b"H[+]" in self.list_components: + self.id_h = self.id_component(b"H[+]") + self.get_ph = self.get_ph_main; + else: + self.id_h = self.id_aqueous(b"H[+]") + self.id_logh = len(self.list_components) + self.id_h + self.get_ph = self.get_ph_second + + def parse_gas_nodes(self): + liq_sat = self.get_liquid_saturation_vs_x(self.ts[0]) + self.is_gas_node = [(not s > 0) for s in liq_sat] + + + def get_timestep(self, timestep): + return self.f[self.map_ts[timestep]] + + def id_species(self, species_list, label): + """ Return the id of a species.""" + return species_list.index(label) + + def id_component(self, label): + """ Return the id of a component.""" + return self.id_species(self.list_components, label) + + def id_mineral(self, label): + """ Return the id of a mineral.""" + return self.id_species(self.list_minerals, label) + + def id_mineral_in_main(self, label): + """ Return the id of a mineral in the main variables vector.""" + return len(self.list_components)+1+self.id_mineral(label) + + def id_gas(self, label): + """ Return the id of a gas.""" + return self.id_species(self.list_gas, label) + + def id_aqueous(self, label): + """ Return the id of a secondary aqueous species.""" + return self.id_species(self.list_aqueous, label) + + # Diffusion profile + # ================= + + def get_diffusion_profile_fit(self, ydata): + """Fit a diffusion profile""" + sqrt_t = np.sqrt(self.ts) + (popt, pcov) = sco.curve_fit(diffusion_profile, sqrt_t, ydata) + return (popt[0], np.sqrt(pcov[0,0])) + + # Carbonation depth + # ================= + + def find_last_true(self, pred): + """Find the last node where pred is true""" + current = 0 + for node in self.nodes: + if pred(node): + current=node + else: + break + return current + + def get_carbonation_depth(self, timestep, pH_threshold=9): + """Return the carbonation depth""" + node = self.find_last_true(lambda n: self.get_ph(timestep, n) < 9) + return self.coordinates[node] + + def get_carbonation_depth_vs_t(self, pH_threshold=9): + """Return the carbonation depths""" + return [self.get_carbonation_depth(t, pH_threshold) for t in self.ts] + + # Main variables + # ============== + + def get_main_variable_vs_x(self, timestep, component, variable): + """ Return the values of a main variable at a timestep.""" + return self.get_timestep(timestep)['main_variables'][component][variable][:] + + def get_main_variable_vs_t(self, node, component, variable): + """ Return the values of a main variable at a node.""" + return [self.f[t]['main_variables'][component][variable][node] for t in self.str_ts] + + # Transport porosity + # ================== + + def get_transport_property_vs_x(self, timestep, property): + """ Return a transport property at a timestep. """ + return self.get_timestep(timestep)['transport_properties'][property][:] + + def get_transport_property_vs_t(self, node, property): + """ Return the values of a transport property at a node. """ + return [self.f[t]['transport_properties'][property][node] for t in self.str_ts] + + # Solid phase volume fraction + # =========================== + + def get_mineral_volfrac(self, timestep, node, mineral): + """Return the volume fraction of a mineral""" + if self.is_gas_node[node]: + return 0 + return self.get_timestep(timestep)[ + 'chemistry_solution'][str(node)][ + 'main_variables'][self.id_mineral_in_main(mineral)] + + + def get_mineral_volfrac_vs_x(self, timestep, mineral): + """ Return the values of the volume fraction of a mineral at a timestep. """ + return [self.get_mineral_volfrac(timestep, node, mineral) + for node in self.nodes] + + def get_mineral_volfrac_vs_t(self, node, mineral): + """ Return the values of the volume fraction of a mineral at a node. """ + return [self.get_mineral_volfrac(timestep, node, mineral) + for timestep in self.ts] + + + # Liquid saturation + # ================= + + def get_liquid_saturation_vs_x(self, timestep): + """ Return the values of the liquid saturation.""" + return self.get_main_variable_vs_x(timestep, b'H2O', 'liquid_saturation') + + def get_liquid_saturation_vs_t(self, node): + """ Return the values of the liquid saturation at node. """ + return self.get_main_variable_vs_t(node, b'H2O', 'liquid_saturation') + + def get_porosity(self, timestep): + """ Return the values of the porosity. """ + return self.get_transport_properties(timestep, 'porosity') + + def get_profile_porosity(self, node): + """ Return the values of the porosity at node. """ + return self.get_profile_transport_properties(node, 'porosity') + + def plot_profile_porosity(self, nodes, **kargs): + """ Plot the profiles of porosity""" + if hasattr(nodes, "__iter__"): + for node in nodes: + plt.plot(self.xs, self.get_profile_porosity(node),label="Node : {0}".format(node)) + else: + plt.plot(self.xs, self.get_profile_porosity(nodes),label="Node : {0}".format(node)) + plt.xlabel("Time (s)") + plt.ylabel("Porosity") + plt.legend() + plt.show() + + def get_ph_vs_t(self, node): + """Return the pH as a function of time""" + return [self.get_ph(t, node) for t in self.ts] + + def get_ph_vs_x(self, timestep): + """Return the pH as a function of space""" + return [self.get_ph(timestep, node) for node in self.nodes] + + + def get_diff(self, t1, t2, func, *args): + """Return the difference between two variables""" + return np.array(func(t2, *args)) - np.array(func(t1, *args)) + +class UnsaturatedPlotter: + """Plot graphs""" + def __init__(self, file): + self.file = file + ###FIXME + self.length_unit="cm" + self.conc_unit="mol/cm^3" + self.show = True + + def plot_oneormore(self, axis, xs, ys_getter, param, label_getter): + """Helper function to plot one or more curve on the same plot""" + if hasattr(param, "__iter__"): + for p in param: + axis.plot(xs, ys_getter(p), label=label_getter(p)) + else: + axis.plot(xs, ys_getter(param), label=label_getter(param)) + + def plot_multiple(self, params, plot_f, *args, **kwargs): + """Plot similar plots in subplots""" + nb_plots = len(params) + fig, axis = plt.subplots(nb_plots) + show_save = self.show; + self.show = False + for it, p in enumerate(params): + plot_f(p, *args, axis=axis[it], **kwargs) + self.show = show_save + if self.show: + plt.show() + + # Solid phases + # ============ + + def plot_mineral_vs_t(self, node, minerals, axis=None): + """Plot solid phases vs t for a node""" + f = self.file + if axis is None: + fig, ax = plt.subplots() + else: + ax = axis + self.plot_oneormore( + ax, f.ts, + lambda m: f.get_mineral_volfrac_vs_t(node, m), + minerals, + lambda m: m + ) + ax.set_xlabel("Time (s)") + ax.set_ylabel("Volume Fraction") + ax.set_title("Node {0}".format(node)) + ax.legend() + plt.legend() + if self.show: + plt.show() + + def plot_mineral_vs_x(self, timestep, minerals, axis=None): + """Plot solid phases vs x for a timestep""" + f = self.file + if axis is None: + fig, ax = plt.subplots() + else: + ax = axis + self.plot_oneormore( + ax, f.coordinates, + lambda m: f.get_mineral_volfrac_vs_x(timestep, m), + minerals, + lambda m: m + ) + ax.set_xlabel("x ({0})".format(self.length_unit)) + ax.set_ylabel("Volume Fraction") + ax.set_title("T = {0} h".format(timestep/3600)) + ax.legend() + if self.show: + plt.show() + + # Main variable + # ============= + + def plot_main_variable_vs_x(self, timestep, component, variable + , axis=None + , ylabel="" + ): + """ Plot the profile for different components at a node """ + f = self.file + if axis is None: + fig, ax = plt.subplots() + else: + ax = axis + if not hasattr(timestep, "__iter__"): + self.plot_oneormore( + ax, f.coordinates, + lambda c: f.get_main_variable_vs_x(timestep, c, variable), + component, + lambda c: c + ) + ax.set_xlabel("x ({0})".format(self.length_unit)) + ax.set_ylabel(ylabel) + ax.set_title("T = {0} h".format(timestep/3600)) + else: + self.plot_oneormore( + ax, f.coordinates, + lambda t: f.get_main_variable_vs_x(t, component, variable), + timestep, + lambda t: "{0} h".format(t/3600) + ) + ax.set_xlabel("x ({0})".format(self.length_unit)) + ax.set_ylabel(ylabel) + ax.set_title("{0} - {1}".format(component, variable)) + ax.legend() + if self.show: + plt.show() + + + def plot_main_variable_vs_t(self, node, component, variable + , axis=None + , ylabel="" + ): + """ Plot the profile for different components at a node """ + f = self.file + if axis is None: + fig, ax = plt.subplots() + else: + ax = axis + if not hasattr(node, "__iter__"): + self.plot_oneormore( + ax, f.ts, + lambda c: f.get_main_variable_vs_t(node, c, variable), + component, + lambda c: c + ) + ax.set_xlabel("t (s)") + ax.set_ylabel(ylabel) + ax.set_title("x={0} {1}".format(f.coordinates[node],self.length_unit)) + ax.legend() + else: + self.plot_oneormore( + ax, f.ts, + lambda n: f.get_main_variable_vs_t(n, component, variable), + node, + lambda n: f.coordinates[n] + ) + ax.set_xlabel("t (s)") + ax.set_ylabel(ylabel) + ax.set_title(component) + ax.legend() + if self.show: + plt.show() + + # liquid saturation + # ----------------- + + def plot_liquid_saturation_vs_x(self, timestep + , axis=None + ): + """ Plot the profile for different components at a node """ + self.plot_main_variable_vs_x(timestep, b"H2O" + , "liquid_saturation" + , axis=axis + , ylabel="liquid saturation" + ) + + + def plot_liquid_saturation_vs_t(self, node + , axis=None + ): + """Plot the profile for different components at a node """ + self.plot_main_variable_vs_t(node, b"H2O" + , "liquid_saturation" + , axis=axis + , ylabel="liquid saturation" + ) + + # Capillary pressure + # ------------------ + + def plot_capillary_pressure_vs_x(self, timestep + , axis=None + ): + """ Plot the profile for different components at a node """ + self.plot_main_variable_vs_x(timestep, b"H2O" + , "capillary_pressure" + , axis=axis + , ylabel="Capillary Pressure (Pa)" + ) + + + def plot_capillary_pressure_vs_t(self, node + , axis=None + ): + """ Plot the profile for different components at a node """ + self.plot_main_variable_vs_t(node, b"H2O" + , "capillary_pressure" + , axis=axis + , ylabel="Capillary Pressure (Pa)" + ) + + # Solid concentration + # ------------------- + + def plot_solid_concentration_vs_x(self, timestep, components + , axis=None + ): + """ Plot the profile for different components at a node """ + self.plot_main_variable_vs_x(timestep, components + , "solid_concentration" + , axis=axis + , ylabel="Solid concentration ({0})".format(self.conc_unit) + ) + + + def plot_solid_concentration_vs_t(self, node, components + , axis=None + ): + """ Plot the profile for different components at a node """ + self.plot_main_variable_vs_t(node, components + , "solid_concentration" + , axis=axis + , ylabel="Solid concentration ({0})".format(self.conc_unit) + ) + + # Aqueous concentration + # ---------------------- + + def plot_aqueous_concentration_vs_x(self, timestep, components + , axis=None + ): + """ Plot the profile for different components at a node """ + self.plot_main_variable_vs_x(timestep, components + , "aqueous_concentration" + , axis=axis + , ylabel="Aqueous concentration ({0})".format(self.conc_unit) + ) + + + def plot_aqueous_concentration_vs_t(self, node, components + , axis=None + ): + """ Plot the profile for different components at a node """ + self.plot_main_variable_vs_t(node, components + , "aqueous_concentration" + , axis=axis + , ylabel="Aqueous concentration ({0})".format(self.conc_unit) + ) + + # Partial pressure + # ---------------- + + def plot_pressure_vs_x(self, timestep, components + , axis=None + ): + """ Plot the profile for different components at a node """ + self.plot_main_variable_vs_x(timestep, components + , "partial_pressure" + , axis=axis + , ylabel="Partial pressure (Pa)" + ) + + def plot_pressure_vs_t(self, node, components + , axis=None + ): + """ Plot the profile for different components at a node """ + self.plot_main_variable_vs_t(node, components + , "partial_pressure" + , axis=axis + , ylabel="Partial pressure (Pa)" + ) + + # pH + # === + + def plot_ph_vs_x(self, timesteps, axis=None): + """Plot the pH as function of x""" + f = self.file + if axis is None: + fig, ax = plt.subplots() + else: + ax = axis + self.plot_oneormore( + ax, f.coordinates, + lambda t: f.get_ph_vs_x(t), + timesteps, + lambda t: "T = {0} h".format(t/3600) + ) + ax.set_xlabel("x ({0})".format(self.length_unit)) + ax.set_ylabel("pH") + ax.set_title("pH") + ax.legend() + if self.show: + plt.show() + + + def plot_ph_vs_t(self, nodes, axis=None): + """Plot the pH as a function of t""" + f = self.file + if axis is None: + fig, ax = plt.subplots() + else: + ax = axis + self.plot_oneormore( + ax, f.ts, + lambda n: f.get_ph_vs_t(n), + nodes, + lambda n: "{0}".format(f.coordinates[n]) + ) + ax.set_xlabel("Time (s)") + ax.set_ylabel("pH") + ax.set_title("pH") + ax.legend() + if self.show: + plt.show() + + # transport property + # ================== + + + def plot_transport_property_vs_x(self, timestep, variable + , axis = None + , ylabel = "" + ): + """Plot the transport property as a function of x""" + f = self.file + if axis is None: + fig, ax = plt.subplots() + else: + ax = axis + if not hasattr(timestep, "__iter__"): + self.plot_oneormore( + ax, f.coordinates, + lambda v: f.get_transport_property_vs_x(timestep, v), + variable, + lambda v: v + ) + ax.set_xlabel("x ({0})".format(self.length_unit)) + ax.set_ylabel(ylabel) + ax.set_title("T = {0} h".format(timestep/3600)) + else: + self.plot_oneormore( + ax, f.coordinates, + lambda t: f.get_transport_property_vs_x(t, variable), + timestep, + lambda t: "{0} h".format(t/3600) + ) + ax.set_xlabel("x ({0})".format(self.length_unit)) + ax.set_ylabel(ylabel) + ax.set_title("{0}".format(variable)) + ax.legend() + if self.show: + plt.show() + + def plot_transport_property_vs_t(self, node, variable + , axis = None + , ylabel = None + ): + """Plot a transport property as a function of t""" + f = self.file + if axis is None: + fig, ax = plt.subplots() + else: + ax = axis + if not hasattr(node, "__iter__"): + self.plot_oneormore( + ax, f.ts, + lambda v: f.get_transport_property_vs_t(node, v), + variable, + lambda v: v + ) + ax.set_xlabel("t (s)") + ax.set_ylabel(ylabel) + ax.set_title("x={0} {1}".format(f.coordinates[node],self.length_unit)) + ax.legend() + else: + self.plot_oneormore( + ax, f.ts, + lambda n: f.get_transport_property_vs_t(n, variable), + node, + lambda n: f.coordinates[n] + ) + ax.set_xlabel("t (s)") + ax.set_ylabel(ylabel) + ax.set_title(variable) + ax.legend() + if self.show: + plt.show() + + def plot_porosity_vs_x(self, timestep): + """Plot the porosity as a function of x""" + self.plot_transport_property_vs_x( + timestep, "porosity", ylabel="porosity", axis=None) + + def plot_porosity_vs_t(self, node): + """Plot the porosity as a function of t""" + self.plot_transport_property_vs_t( + node, "porosity", ylabel="porosity", axis=None) + + def plot_carbonation_depth_vs_t(self, pH_threshold=9): + """Plot the carbonation depth as a function of t""" + f = self.file + plt.plot(f.ts, f.get_carbonation_depth_vs_t(pH_threshold), "o") + plt.xlabel("Time (s)".format(self.length_unit)) + plt.ylabel("Carbonation depth ({0})".format(self.length_unit)) + if self.show: + plt.show() + + def plot_diffusion_fit(self, ydata, ylabel=""): + """Plot a diffusion profile fit""" + f = self.file + sqrt_t = np.sqrt(f.ts) + (a, stda) = f.get_diffusion_profile_fit(ydata); + fit = [diffusion_profile(sq, a) for sq in sqrt_t] + plt.plot(sqrt_t, ydata, "+", label="simulation") + plt.plot(sqrt_t, fit, label="a={0} +/- {1}".format(a, stda)) + plt.xlabel("sqrt(T) (sqrt(s))") + plt.ylabel(ylabel) + plt.legend() + if self.show: + plt.show() + # open the file print(" * load file") -try: - filename = sys.argv[1] - f = h5py.File(filename, mode='r') -except Exception as e: - fail("Error while opening file : {0}".format(filename), e) - print(" * analyse file") -# timestep -# ======== - -try: - str_xs = [x for x in f][:-2] - str_xs.sort(key=float) - xs = [float(x) for x in str_xs] - map_xs = dict(zip(xs, str_xs)) -except Exception as e: - fail("Error while parsing timesteps", e) - -# database -# ========= - -try: - database = f['database'] -except Exception as e: - fail("Error while parsing database section", e) - -def get_list_species(species): - return list(database[species][:]) - -try: - list_components = get_list_species('components') - list_minerals = get_list_species('minerals') - list_gas = get_list_species('gas') - list_aqueous = get_list_species('aqueous') -except Exception as e: - fail("Error while parsing database labels", e) - - -def get_timestep(timestep): - return f[map_xs[timestep]] - -def id_species(species_list, label): - """ Return the id of a species.""" - return species_list.index(label) - -def id_component(label): - """ Return the id of a component.""" - return id_species(list_components, label) - -def id_mineral(label): - """ Return the id of a mineral.""" - return id_species(list_minerals, label) - -def id_mineral_in_main(label): - """ Return the id of a mineral in the main variables vector.""" - return len(list_components)+1+id_mineral(label) - -def id_gas(label): - """ Return the id of a gas.""" - return id_species(list_gas, label) - -def id_aqueous(label): - """ Return the id of a secondary aqueous species.""" - return id_species(list_aqueous, label) - - - -# mesh -# ==== - -try: - mesh = f['mesh'] - coordinates = mesh['coordinates'][:] - nb_nodes = coordinates.shape[0] - nodes = (int(x) for x in range(nb_nodes)) -except Exception as e: - fail("Error while parsing mesh", e) - -# variables -# ========= - -def get_main_variables(timestep, component, variable): - """ Return the values of a main variable at a timestep.""" - return get_timestep(timestep)['main_variables'][component][variable][:] - -def get_profile_main_variables(node, component, variable): - """ Return the values of a main variable at a node.""" - return [f[x]['main_variables'][component][variable][node] for x in str_xs] - -def plot_profile_components(node, components, variable): - """ Plot the profile for different components at a node """ - if hasattr(components, "__iter__"): - for component in components: - plt.plot(xs, - get_profile_main_variables(node, component, variable), - label="{0}".format(component) - ) - else: - plt.plot(xs, - get_profile_main_variables(node, components, variable), - label="{0}".format(components) - ) - plt.xlabel("Time (s)") - plt.ylabel("{0}, node {1}".format(variable, node)) - plt.legend() - plt.show() - -def plot_profile_component_nodes(nodes, component, variable): - """ Plot the profile for a component at nodes """ - if hasattr(nodes, "__iter__"): - for node in nodes: - plt.plot(xs, - get_profile_main_variables(node, component, variable), - label="Node {0}".format(node) - ) - else: - plt.plot(xs, - get_profile_main_variables(nodes, components, variable), - label="Node {0}".format(node) - ) - plt.xlabel("Time (s)") - plt.ylabel("{0} {1}".format(variable, component)) - plt.legend() - plt.show() - -def get_transport_properties(timestep, prop): - """ Return a transport properties at a timestep. """ - return get_timestep(timestep)['transport_properties'][prop][:] - -def get_profile_transport_properties(node, prop): - """ Return the values of a transport properties at a node. """ - return [f[x]['transport_properties'][prop][node] for x in str_xs] - -def get_mineral_volfrac(timestep, mineral): - """ Return the values of the volume fraction of a mineral at a timestep. """ - return [get_timestep(timestep)[ - 'chemistry_solution'][str(node)][ - 'main_variables'][id_mineral_in_main(mineral)] - for node in nodes] - -def get_profile_mineral_volfrac(node, mineral): - """ Return the values of the volume fraction of a mineral at a node. """ - return [f[x][ - 'chemistry_solution'][str(node)][ - 'main_variables'][id_mineral_in_main(mineral)] - for x in str_xs] - -def get_liquid_saturation(timestep): - """ Return the values of the liquid saturation.""" - return get_main_variables(timestep, b'H2O', 'liquid_saturation') - -def get_profile_liquid_saturation(node): - """ Return the values of the liquid saturation at node. """ - return get_profile_main_variables(node, b'H2O', 'liquid_saturation') - -def plot_profile_liquid_saturation(nodes, **kargs): - """ Plot the profile of liquid saturation """ - if hasattr(nodes, "__iter__"): - for node in nodes: - plt.plot(xs, get_profile_liquid_saturation(node),label="Node : {0}".format(node)) - else: - plt.plot(xs, get_profile_liquid_saturation(nodes),label="Node : {0}".format(nodes)) - plt.xlabel("Time (s)") - plt.ylabel("Liquid saturation") - plt.legend() - plt.show() - -def get_porosity(timestep): - """ Return the values of the porosity. """ - return get_transport_properties(timestep, 'porosity') - -def get_profile_porosity(node): - """ Return the values of the porosity at node. """ - return get_profile_transport_properties(node, 'porosity') - -def plot_profile_porosity(nodes, **kargs): - """ Plot the profiles of porosity""" - if hasattr(nodes, "__iter__"): - for node in nodes: - plt.plot(xs, get_profile_porosity(node),label="Node : {0}".format(node)) - else: - plt.plot(xs, get_profile_porosity(nodes),label="Node : {0}".format(node)) - plt.xlabel("Time (s)") - plt.ylabel("Porosity") - plt.legend() - plt.show() - -def plot_mineral_profile(node): - for mineral in list_minerals: - plt.plot(xs, get_profile_mineral_volfrac(node, mineral), label=mineral); - plt.legend() - plt.show() - -def plot_mineral_profiles(nodes, mineral): - if hasattr(nodes, "__iter__"): - for node in nodes: - plt.plot(xs, get_profile_mineral_volfrac(node, mineral), label="Node {0}".format(node)) - else: - plt.plot(xs, get_profile_mineral_volfrac(node, mineral), label="Node {0}".format(node)) - plt.xlabel("Time (s)") - plt.ylabel("Volume fraction {0}".format(mineral)) - plt.legend() - plt.show() - - -if b"H[+]" in list_components: - id_h = id_component(b"H[+]") - id_logh = id_h - def get_ph(timestep, node): - chem_sol = get_timestep(timestep)['chemistry_solution'][str(node)] - log_h = chem_sol['main_variables'][id_h] - log_gamma_h = chem_sol['log_gamma'][id_h] - return -(log_gamma_h + log_h) -else: - id_h = id_aqueous(b"H[+]") - id_logh = len(list_components) + id_h - def get_ph(timestep, node): - chem_sol = get_timestep(timestep)['chemistry_solution'][str(node)] - log_h = np.log10(chem_sol['secondary_molalities'][id_h]) - log_gamma_h = chem_sol['log_gamma'][id_h] - return -(log_gamma_h + log_h) - -def get_profile_ph(node): - return [get_ph(x, node) for x in xs] - - -def plot_ph_profiles(nodes): - if hasattr(nodes, "__iter__"): - for node in nodes: - plt.plot(xs, get_profile_ph(node), label="Node {0}".format(node)) - else: - plt.plot(xs, get_profile_ph(node), label="Node {0}".format(node)) - plt.xlabel("Time (s)") - plt.ylabel("pH") - plt.legend() - plt.show() +sol = UnsaturatedHDF5File(sys.argv[1]); +plotter = UnsaturatedPlotter(sol) print(" * starting interactive session")