Page MenuHomec4science
Files F92216702

analyse_unsaturated_hdf5.py
No OneTemporary
Actions

File Metadata

Created
Mon, Nov 18, 10:25

analyse_unsaturated_hdf5.py
View Options

# -----------------------------------------------------------------------------
#Copyright (c) 2015 Fabien Georget <fabieng@princeton.edu>, 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[+]")
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_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")

Event Timeline

c4science · Help