Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/scripts/analyser/analyse_unsaturated_hdf5.py
06af08aa99cemaster

Last Change
Raw File

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
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")
print(" * analyse file")
sol = UnsaturatedHDF5File(sys.argv[1]);
plotter = UnsaturatedPlotter(sol)
print(" * starting interactive session")
c4science · Help