diff --git a/scripts/CMakeLists.txt b/scripts/CMakeLists.txt index e451361..097cd3f 100644 --- a/scripts/CMakeLists.txt +++ b/scripts/CMakeLists.txt @@ -1,35 +1,52 @@ add_custom_target(scripts_incl SOURCES FindSpecMiCP.cmake reactmicp.qsub pgo_sequence.sh specmicp.pc.in reactmicp.pc.in + + analyser/reactmicp_unsaturated_analyse.sh + analyser/analyse_unsaturated_hdf5.py ) # the cmake find module # --------------------- install(FILES FindSpecMiCP.cmake DESTINATION ${SHARE_INSTALL_DIR}/cmake ) # openPBS submission file file(COPY reactmicp.qsub DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) install(FILES reactmicp.qsub DESTINATION ${SHARE_INSTALL_DIR} ) # Profile guided optimization file(COPY pgo_sequence.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) # pkg-config files configure_file(specmicp.pc.in ${CMAKE_CURRENT_BINARY_DIR}/specmicp.pc @ONLY) install(FILES ${CMAKE_CURRENT_BINARY_DIR}/specmicp.pc DESTINATION ${LIBRARY_INSTALL_DIR}/pkgconfig ) configure_file(reactmicp.pc.in ${CMAKE_CURRENT_BINARY_DIR}/reactmicp.pc @ONLY) install(FILES ${CMAKE_CURRENT_BINARY_DIR}/reactmicp.pc DESTINATION ${LIBRARY_INSTALL_DIR}/pkgconfig ) + +# analyser + +file(COPY + analyser/reactmicp_unsaturated_analyse.sh + analyser/analyse_unsaturated_hdf5.py + DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/analyser +) + +install(FILES + analyser/reactmicp_unsaturated_analyse.sh + analyser/analyse_unsaturated_hdf5.py + DESTINATION ${SHARE_INSTALL_DIR}/analyser +) diff --git a/scripts/analyser/analyse_unsaturated_hdf5.py b/scripts/analyser/analyse_unsaturated_hdf5.py new file mode 100644 index 0000000..4da1d4d --- /dev/null +++ b/scripts/analyser/analyse_unsaturated_hdf5.py @@ -0,0 +1,182 @@ +# ----------------------------------------------------------------------------- +#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 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 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_mineral_profile(node): + for mineral in list_minerals: + plt.plot(xs, get_profile_mineral_volfrac(node, mineral), label=mineral); + plt.legend() + plt.show() + + + +print(" * starting interactive session") diff --git a/scripts/analyser/reactmicp_unsaturated_analyse.sh b/scripts/analyser/reactmicp_unsaturated_analyse.sh new file mode 100755 index 0000000..ca1f09b --- /dev/null +++ b/scripts/analyser/reactmicp_unsaturated_analyse.sh @@ -0,0 +1,194 @@ +#! /usr/bin/env bash + +# ----------------------------------------------------------------------------- +# 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. +# ----------------------------------------------------------------------------- + +# This file is a wrapper script to analyse a hdf5 file produced by the +# unsaturated system of ReactMiCP +# It calls a python interpreter in interactive mode and loads a module + + +# The module where the useful functions are defined +python_script=analyse_unsaturated_hdf5.py + +# Return codes +RETCODE_NO_FILE=1 +RETCODE_NO_PYTHON=2 +RETCODE_FILE_DOESNT_EXIST=3 +RETCODE_BAD_OPTION=4 +RETCODE_BAD_PYTHON_MODULE=5 + +echo " +=================================== + + ReactMiCP + =========== + + Unsaturated HDF5 output analyser + -------------------------------- + +=================================== + +Copyright (c) 2015 F. Georget + +" + +function check_file_arg { + # check that the file argument exist + if [[ ! $1 || -z $1 ]] + then + echo "Error : one argument required. Aborting" >&2; + exit ${RETCODE_NO_FILE}; + fi +} + +function print_help { + echo " +Analyse hdf5 file produced by the unsaturated system of ReactMiCP. +The analyser is a script which open a python interpreter and load useful functions. + +Usage : $0 [options] + + : hdf5 file + +Options : + -m, --module alternative python module to use + -p, --python python interpreter to use + -h, --h : print this help" + exit 0 +} + +check_file_arg $1 + +# Options +# ------- +while [[ $# > 1 ]] +do + key="$1" + case ${key} in + -h|--help) + print_help + ;; + -m|--module) + python_script="$2" + shift + ;; + -p|--python) + python_shell="$2" + shift + ;; + *) + echo "Unknow option '$1'. Aborting" >&2 + exit ${RETCODE_BAD_OPTION} + ;; + esac +shift +done + +check_file_arg $1 + +if [[ $1 = "-h" || $1 = "--help" ]] +then + print_help +fi + +# check that the HDF5 file exist +# no point going further if it doesn't +if [[ ! -f $1 ]] +then + echo "Error : no such file : '$1'. Aborting" >&2; + exit ${RETCODE_FILE_DOESNT_EXIST} +fi + +# copy file parameter +filepath=$1 +shift + +# Find the python script +if [[ ! -f ${python_script} ]] +then + # Find the true directory of this bash script + # https://stackoverflow.com/questions/59895/can-a-bash-script-tell-what-directory-its-stored-in + SOURCE="${BASH_SOURCE[0]}" + while [ -h "$SOURCE" ]; do # resolve $SOURCE until the file is no longer a symlink + DIR="$( cd -P "$( dirname "$SOURCE" )" && pwd )" + SOURCE="$(readlink "$SOURCE")" + [[ $SOURCE != /* ]] && SOURCE="$DIR/$SOURCE" # if $SOURCE was a relative symlink, we need to resolve it relative to the path where the symlink file was located + done + DIR="$( cd -P "$( dirname "$SOURCE" )" && pwd )" + + test_src=${DIR}/${python_script} + if [[ ! -f ${test_src} ]] + then + echo "Unable to find python module. Aborting" >&2; + exit ${RETCODE_BAD_PYTHON_MODULE}; + else + python_script=${test_src} + fi +fi + + +# find the python shell +if [[ -z ${python_shell} ]] +then + python_shell=$(command -v ipython) + if [[ $? ]] + then + echo " * python shell detected : ipython (${python_shell})" + else + python_shell=$(command -v python) + if [[ $1 ]] + then + echo " * python shell detected : python (${python_shell})" + else + echo "Error : A python interpreter is required. Aborting\n" >&2; + exit ${RETCODE_NO_PYTHON}; + fi + fi +fi + +# Don't print the python banner if possible +# it makes thing just more confusing +basename_python=$(basename ${python_shell}) +if [[ ${basename_python} = ipython* ]] +then + no_banner="--no-banner" +elif [[ ${basename_python} = python* ]] +then + python_version=$(${python_shell} -V 2>&1 | cut -d' ' -f2 | cut -d'.' -f1 -) + if [[ $python_version -ge 3 ]] + then + no_banner="-q" + fi +fi + +# run the analysis +exec ${python_shell} -i ${no_banner} ${python_script} ${filepath}