diff --git a/scripts/analyser/analyse_unsaturated_hdf5.py b/scripts/analyser/analyse_unsaturated_hdf5.py index 4da1d4d..be82b7f 100644 --- a/scripts/analyser/analyse_unsaturated_hdf5.py +++ b/scripts/analyser/analyse_unsaturated_hdf5.py @@ -1,182 +1,250 @@ # ----------------------------------------------------------------------------- #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 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() print(" * starting interactive session")