diff --git a/scripts/analyser/analyse_unsaturated_hdf5.py b/scripts/analyser/analyse_unsaturated_hdf5.py index 11c27ed..861c1b9 100644 --- a/scripts/analyser/analyse_unsaturated_hdf5.py +++ b/scripts/analyser/analyse_unsaturated_hdf5.py @@ -1,287 +1,287 @@ # ----------------------------------------------------------------------------- #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 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) # 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 "H[+]" in list_components: - id_h = id_component("H[+]") +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("H[+]") + 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() print(" * starting interactive session")