diff --git a/src/reactmicp/io/hdf5_unsaturated.cpp b/src/reactmicp/io/hdf5_unsaturated.cpp index 45fb569..6e03741 100644 --- a/src/reactmicp/io/hdf5_unsaturated.cpp +++ b/src/reactmicp/io/hdf5_unsaturated.cpp @@ -1,456 +1,746 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. 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. * ============================================================================= */ #include "hdf5_unsaturated.hpp" #include "reactmicp/systems/unsaturated/variables.hpp" -#include -#include "specmicp_common/io/specmicp_hdf5.hpp" -#include "specmicp_common/io/hdf5_eigen.hpp" +#include "specmicp_common/io/hdf5/file.hpp" +#include "specmicp_common/io/hdf5/group.hpp" +#include "specmicp_common/io/hdf5/dataset.hpp" #include "specmicp_common/io/hdf5_timesteps.hpp" + #include "specmicp/io/hdf5_adimensional.hpp" #include "dfpm/io/hdf5_mesh.hpp" #include "specmicp_database/io/hdf5_database.hpp" #include "specmicp_common/physics/units.hpp" #include "specmicp_database/data_container.hpp" #include "specmicp_common/compat.hpp" +#include "specmicp_common/log.hpp" #include #include +#include using namespace specmicp::reactmicp::systems::unsaturated; #define MESH_GRPNAME "mesh" #define DATABASE_GRPNAME "database" #define UNIT_ATTRIBUTE "unit" #define MAIN_VARIABLES_GRPNAME "main_variables" #define CHEMISTRY_SOLUTION_GRPNAME "chemistry_solution" #define TRANSPORT_PROPERTIES_GRPNAME "transport_properties" #define LIQUID_SATURATION_DSET "liquid_saturation" #define AQUEOUS_CONC_DSET "aqueous_concentration" #define SOLID_CONC_DSET "solid_concentration" #define PRESSURE_DSET "partial_pressure" #define CAP_PRESSURE_DSET "capillary_pressure" #define POROSITY_DSET "porosity" #define LIQUID_DIFFUSIVITY_DSET "liquid_diffusivity" #define REL_LIQUID_DIFFUSIVITY_DSET "relative_liquid_diffusivity" #define LIQUID_PERMEABILITY_DSET "liquid_permeability" #define REL_LIQUID_PERMEABILITY_DSET "relative_liquid_permeability" #define RESISTANCE_GAS_DIFFUSIVITY_DSET "resistance_gas_diffusivity" #define REL_GAS_DIFFUSIVITY_DSET "relative_gas_diffusivity" namespace specmicp { namespace io { //! \brief Implementation of the UnsaturatedHDF5Saver //! //! \internal struct SPECMICP_DLL_LOCAL UnsaturatedHDF5Saver::UnsaturatedHDF5SaverImpl { public: std::string m_filepath; std::shared_ptr m_vars; //!< the variables database::RawDatabasePtr m_database; //!< the database UnsaturatedHDF5SaverImpl( const std::string& filename, std::shared_ptr vars, const units::UnitsSet the_units): m_filepath(filename), m_vars(vars), m_database(vars->get_database()) { init_file(the_units); } //! \brief Save a timestep void save_timestep( - const std::string& name, - const std::string& section + const std::string& name ); //! \brief Save the main variables - //! - //! \param file HDF5 file - //! \param section section where to save the main variables void save_main_variables( - io::HDF5File& file, - const std::string& section + hdf5::Group& timestep_grp ); //! \brief Save the chemical solution - //! - //! \param file HDF5 file - //! \param section section where to save the chemical solutions void save_chemical_solution( - io::HDF5File& file, - const std::string& section + hdf5::Group& timestep_grp ); //! \brief Save the transport properties - //! - //! \param file HDF5 file - //! \param section section where to save the transport properties void save_transport_properties( - io::HDF5File& file, - const std::string& section + hdf5::Group& timestep_grp ); + //! \brief create the file and save common values (mesh, db, units) void init_file(const units::UnitsSet& the_units); }; -struct SPECMICP_DLL_LOCAL UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl -{ -public: - UnsaturatedHDF5ReaderImpl(const std::string& name): - m_file(name, HDF5_OpenMode::OpenReadOnly), - m_timesteps(m_file) - { - } - - void initialize_variables( - scalar_t timestep, - reactmicp::systems::unsaturated::UnsaturatedVariables* vars - ); - - void read_main_variables( - const H5::Group& timestep_group, - reactmicp::systems::unsaturated::UnsaturatedVariables* vars - ); - void read_chemical_solutions( - const H5::Group& timestep_group, - reactmicp::systems::unsaturated::UnsaturatedVariables* vars - ); - void read_transport_properties( - const H5::Group& timestep_group, - reactmicp::systems::unsaturated::UnsaturatedVariables* vars - ); - -private: - HDF5File m_file; - HDF5Timesteps m_timesteps; -}; - // UnsaturatedHDF5Saver::UnsaturatedHDF5Saver( const std::string& filename, std::shared_ptr vars, const units::UnitsSet& the_units ): m_impl(make_unique(filename, vars, the_units)) {} UnsaturatedHDF5Saver::~UnsaturatedHDF5Saver() = default; void UnsaturatedHDF5Saver::save_timestep(scalar_t timestep) { - return m_impl->save_timestep(std::to_string(timestep), "/"); + return m_impl->save_timestep(std::to_string(timestep)); } -//! \brief Initialize variables from a saved timestep -//! -//! The Mesh and the database must already be ok -void UnsaturatedHDF5Reader::initialize_variables( - scalar_t timestep, - reactmicp::systems::unsaturated::UnsaturatedVariables* vars - ) -{ - return m_impl->initialize_variables(timestep, vars); -} - // Implementation // ============== //! \brief Save a timestep void UnsaturatedHDF5Saver::UnsaturatedHDF5SaverImpl::save_timestep( - const std::string& name, - const std::string& section + const std::string& name ) { - HDF5File the_file(m_filepath, HDF5_OpenMode::OpenReadWrite); + hdf5::File the_file = hdf5::File::open(m_filepath, hdf5::OpenMode::OpenReadWrite); + hdf5::Group timestep_grp = the_file.create_group(name); - auto group_name = the_file.complete_name(name, section); - auto group = the_file.create_group(group_name); - - save_main_variables(the_file, group_name); - save_chemical_solution(the_file, group_name); - save_transport_properties(the_file, group_name); + save_main_variables(timestep_grp); + save_chemical_solution(timestep_grp); + save_transport_properties(timestep_grp); } +// ugly, but more readable at the end +#define uvars(var_name) \ + m_vars->get_ ## var_name ().variable + +#define uvarsi(var_name, index) \ + m_vars->get_ ## var_name (index).variable + + void UnsaturatedHDF5Saver::UnsaturatedHDF5SaverImpl::save_main_variables( - HDF5File& file, - const std::string& section + hdf5::Group& timestep_grp ) { - auto group_name = file.complete_name(MAIN_VARIABLES_GRPNAME, section); - auto group = file.create_group(group_name); + hdf5::Group main_vars_grp = timestep_grp.create_group(MAIN_VARIABLES_GRPNAME); // water { - auto water_group = group->createGroup(m_database->get_label_component(0)); - - save_eigen_matrix(water_group, LIQUID_SATURATION_DSET, - m_vars->get_liquid_saturation().variable); - - save_eigen_matrix(water_group, AQUEOUS_CONC_DSET, - m_vars->get_water_aqueous_concentration().variable); - - save_eigen_matrix(water_group, SOLID_CONC_DSET, - m_vars->get_solid_concentration(0).variable); - - if (m_vars->component_has_gas(0)) { - save_eigen_matrix(water_group, PRESSURE_DSET, - m_vars->get_pressure_main_variables(0).variable); + auto water_group = main_vars_grp.create_group(m_database->get_label_component(0)); + water_group.create_vector_dataset(LIQUID_SATURATION_DSET, + uvars(liquid_saturation)); + water_group.create_vector_dataset(AQUEOUS_CONC_DSET, + uvars(water_aqueous_concentration)); + water_group.create_vector_dataset(SOLID_CONC_DSET, + uvarsi(solid_concentration, 0)); + if (m_vars->component_has_gas(0)) + { + water_group.create_vector_dataset(PRESSURE_DSET, + uvarsi(pressure_main_variables, 0)); } - - save_eigen_matrix(water_group, CAP_PRESSURE_DSET, - m_vars->get_capillary_pressure().variable); + water_group.create_vector_dataset(CAP_PRESSURE_DSET, + uvars(capillary_pressure)); } for (index_t aq_component: m_database->range_aqueous_component()) { - auto aqcomp_group = group->createGroup(m_database->get_label_component(aq_component)); - - save_eigen_matrix(aqcomp_group, AQUEOUS_CONC_DSET, - m_vars->get_aqueous_concentration(aq_component).variable); - - save_eigen_matrix(aqcomp_group, SOLID_CONC_DSET, - m_vars->get_solid_concentration(aq_component).variable); - - if (m_vars->component_has_gas(aq_component)) { - save_eigen_matrix(aqcomp_group, PRESSURE_DSET, - m_vars->get_pressure_main_variables(aq_component).variable); + auto aqcomp_group = main_vars_grp.create_group(m_database->get_label_component(aq_component)); + + aqcomp_group.create_vector_dataset(AQUEOUS_CONC_DSET, + uvarsi(aqueous_concentration, aq_component)); + aqcomp_group.create_vector_dataset(SOLID_CONC_DSET, + uvarsi(solid_concentration, aq_component)); + if (m_vars->component_has_gas(aq_component)) + { + aqcomp_group.create_vector_dataset(PRESSURE_DSET, + uvarsi(pressure_main_variables, aq_component)); } } } void UnsaturatedHDF5Saver::UnsaturatedHDF5SaverImpl::save_chemical_solution( - HDF5File& file, - const std::string& section + hdf5::Group& timestep_grp ) { - auto group_name = file.complete_name(CHEMISTRY_SOLUTION_GRPNAME, section); - auto group = file.create_group(group_name); - + auto chem_grp = timestep_grp.create_group(CHEMISTRY_SOLUTION_GRPNAME); for (auto node: m_vars->get_mesh()->range_nodes()) { - save_adimensional_system_solution(file, std::to_string(node), group_name, + if (m_vars->get_adim_solution(node).is_valid) + { + save_adimensional_system_solution(chem_grp, + std::to_string(node), m_vars->get_adim_solution(node)); + } } } + + // time to get ugly and use macro // this is simply a quick wrapper that call the save_eigen matrix function // It avoids all the boilerplate #define save_transport_property(name, var) \ - save_eigen_matrix(*group.get(), name, m_vars->get_##var ().variable) + transport_grp.create_vector_dataset(name, uvars(var)) void UnsaturatedHDF5Saver::UnsaturatedHDF5SaverImpl::save_transport_properties( - HDF5File& file, - const std::string& section + hdf5::Group& timestep_grp ) { - auto group_name = file.complete_name(TRANSPORT_PROPERTIES_GRPNAME, section); - auto group = file.create_group(group_name); + auto transport_grp = timestep_grp.create_group(TRANSPORT_PROPERTIES_GRPNAME); save_transport_property(POROSITY_DSET, porosity); save_transport_property(LIQUID_DIFFUSIVITY_DSET, liquid_diffusivity); save_transport_property(REL_LIQUID_DIFFUSIVITY_DSET, relative_liquid_diffusivity); save_transport_property(LIQUID_PERMEABILITY_DSET, liquid_permeability); save_transport_property(REL_LIQUID_PERMEABILITY_DSET, relative_liquid_permeability); save_transport_property(RESISTANCE_GAS_DIFFUSIVITY_DSET, resistance_gas_diffusivity); save_transport_property(REL_GAS_DIFFUSIVITY_DSET, relative_gas_diffusivity); } #undef save_transport_property // we remove the ugly macro, no one saws anything :) - +#undef uvars +#undef uvarsi void UnsaturatedHDF5Saver::UnsaturatedHDF5SaverImpl::init_file( const units::UnitsSet& the_units ) { - HDF5File the_file(m_filepath, HDF5_OpenMode::CreateTruncate); - - save_mesh_coordinates(the_file, MESH_GRPNAME, "/", m_vars->get_mesh()); - save_database_labels(the_file, DATABASE_GRPNAME, "/", m_vars->get_database()); + hdf5::File the_file = hdf5::File::open(m_filepath, hdf5::OpenMode::CreateTruncate); - hsize_t dims[] = {2}; - auto dspace = H5::DataSpace(1, dims); - auto attribute = the_file.create_attribute( - UNIT_ATTRIBUTE, - H5::PredType::NATIVE_DOUBLE, - dspace); + save_mesh_coordinates(the_file, MESH_GRPNAME, m_vars->get_mesh()); + save_database_labels(the_file, DATABASE_GRPNAME, m_vars->get_database()); - scalar_t attribute_values[] = { + std::array units_values = { units::scaling_factor(the_units.length), - units::scaling_factor(the_units.mass) + units::scaling_factor(the_units.mass), + units::scaling_factor(the_units.quantity) }; - attribute->write(H5::PredType::NATIVE_DOUBLE, attribute_values); + the_file.create_scalar_attribute(UNIT_ATTRIBUTE, units_values); } +// =========================================================================== +// =========================================================================== + // Reader // ======= + +struct SPECMICP_DLL_LOCAL UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl +{ + UnsaturatedHDF5ReaderImpl(const std::string& name): + m_file(hdf5::File::open(name, hdf5::OpenMode::OpenReadOnly)), + m_timesteps(m_file) + { + } + + units::UnitsSet get_units(); + + AdimensionalSystemSolution get_adim_solution( + scalar_t timestep, + index_t node + ); + AdimensionalSystemSolution get_adim_solution( + std::string timestep, + index_t node + ); + + void initialize_variables( + scalar_t timestep, + reactmicp::systems::unsaturated::UnsaturatedVariables* vars + ); + + void read_main_variables( + const hdf5::GroupPath& timestep_group, + reactmicp::systems::unsaturated::UnsaturatedVariables* vars + ); + void read_chemical_solutions( + const hdf5::GroupPath& timestep_group, + reactmicp::systems::unsaturated::UnsaturatedVariables* vars + ); + void read_transport_properties( + const hdf5::GroupPath& timestep_group, + reactmicp::systems::unsaturated::UnsaturatedVariables* vars + ); + + //! \brief Return a main variable against x + Vector main_variable_vs_x( + const std::string& timestep, + const std::string& component, + const std::string& variable + ); + //! \brief Return a main variable against t + Vector main_variable_vs_t( + index_t node, + const std::string& component, + const std::string& variable + ); + + //! \brief Open a timestep group + hdf5::Group open_timestep(const std::string& timestep); + //! \brief Open a timestep group + hdf5::Group open_timestep(const scalar_t& timestep); + + //! \brief Open the main variables group + hdf5::Group open_main_variables(const std::string& timestep); + //! \brief Open the main variables group for a component + hdf5::Group open_main_variables_component( + const std::string& timestep, + const std::string& component + ); + + + //! \brief Open the transport variables group + hdf5::Group open_transport_variables(const std::string& timestep); + //! \brief Return a main variable against x + Vector transport_variable_vs_x( + const std::string& timestep, + const std::string& variable + ); + //! \brief Return a main variable against t + Vector transport_variable_vs_t( + index_t node, + const std::string& variable + ); + + hdf5::Group open_chemistry_solutions(const std::string& timestep); + + hdf5::File m_file; + HDF5Timesteps m_timesteps; +}; + + +UnsaturatedHDF5Reader::UnsaturatedHDF5Reader( + const std::string &filepath + ): + m_impl(utils::make_pimpl(filepath)) +{ +} + +UnsaturatedHDF5Reader::~UnsaturatedHDF5Reader() = default; + +//! \brief Initialize variables from a saved timestep +//! +//! The Mesh and the database must already be ok +void UnsaturatedHDF5Reader::initialize_variables( + scalar_t timestep, + reactmicp::systems::unsaturated::UnsaturatedVariables* vars + ) +{ + return m_impl->initialize_variables(timestep, vars); +} + + +const HDF5Timesteps& UnsaturatedHDF5Reader::get_timesteps() const +{ + return m_impl->m_timesteps; +} + +units::UnitsSet UnsaturatedHDF5Reader::get_units() +{ + return m_impl->get_units(); +} + +AdimensionalSystemSolution UnsaturatedHDF5Reader::get_adim_solution( + scalar_t timestep, + index_t node + ) +{ + return m_impl->get_adim_solution(timestep, node); +} +AdimensionalSystemSolution UnsaturatedHDF5Reader::get_adim_solution( + std::string timestep, + index_t node + ) +{ + return m_impl->get_adim_solution(timestep, node); +} + +Vector UnsaturatedHDF5Reader::main_variable_vs_x( + const std::string& timestep, + const std::string& component, + const std::string& variable + ) +{ + return m_impl->main_variable_vs_x(timestep, component, variable); +} +Vector UnsaturatedHDF5Reader::main_variable_vs_x( + scalar_t timestep, + const std::string& component, + const std::string& variable + ) +{ + const std::string str_timestep = m_impl->m_timesteps.get_string(timestep); + return m_impl->main_variable_vs_x(str_timestep, component, variable); +} +Vector UnsaturatedHDF5Reader::main_variable_vs_t( + index_t node, + const std::string& component, + const std::string& variable + ) +{ + return m_impl->main_variable_vs_t(node, component, variable); +} + + +Vector UnsaturatedHDF5Reader::transport_variable_vs_x( + const std::string& timestep, + const std::string& variable + ) +{ + return m_impl->transport_variable_vs_x(timestep, variable); +} +Vector UnsaturatedHDF5Reader::transport_variable_vs_x( + scalar_t timestep, + const std::string& variable + ) +{ + const std::string str_timestep = m_impl->m_timesteps.get_string(timestep); + return m_impl->transport_variable_vs_x(str_timestep, variable); +} + +Vector UnsaturatedHDF5Reader::transport_variable_vs_t( + index_t node, + const std::string& variable + ) +{ + return m_impl->transport_variable_vs_t(node, variable); +} + +// =========== Implementation =============================== + + +units::UnitsSet UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::get_units() +{ + units::UnitsSet the_units; + // quantity conversion factor was not solved in early version + // really early, can probably be simplified at some point + std::vector conversion_factor = m_file.read_scalar_attribute(UNIT_ATTRIBUTE); + the_units.length = units::from_scaling_factor(conversion_factor[0]); + the_units.mass = units::from_scaling_factor(conversion_factor[1]); + if (conversion_factor.size() == 2) + { + WARNING << "Units for quantity of matter was not recorded, assuming moles."; + } + else + { + the_units.quantity = units::from_scaling_factor(conversion_factor[2]); + } + return the_units; +} + + void UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::initialize_variables( scalar_t timestep, reactmicp::systems::unsaturated::UnsaturatedVariables *vars ) { auto name_group = m_timesteps.get_string(timestep); - auto group = m_file.open_group(name_group); + auto timestep_group = m_file.open_group(name_group); - const auto& grp_ref = *group.get(); - read_main_variables(grp_ref, vars); - read_chemical_solutions(grp_ref, vars); - read_transport_properties(grp_ref, vars); + read_main_variables(timestep_group, vars); + read_chemical_solutions(timestep_group, vars); + read_transport_properties(timestep_group, vars); } + +// ugly, but more readable at the end +#define uvars(var_name) \ + vars->get_ ## var_name ().variable + +#define uvarsi(var_name, index) \ + vars->get_ ## var_name (index).variable + + void UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::read_main_variables( - const H5::Group& timestep_group, + const hdf5::GroupPath& timestep_group, reactmicp::systems::unsaturated::UnsaturatedVariables* vars ) { - H5::Group main_vars_grp = timestep_group.openGroup(MAIN_VARIABLES_GRPNAME); + hdf5::Group main_vars_grp = timestep_group.open_group(MAIN_VARIABLES_GRPNAME); database::RawDatabasePtr raw_database = vars->get_database(); { // Water - auto group = main_vars_grp.openGroup(raw_database->get_label_component(0)); - read_eigen_matrix(group, LIQUID_SATURATION_DSET, - vars->get_liquid_saturation().variable); - read_eigen_matrix(group, AQUEOUS_CONC_DSET, - vars->get_water_aqueous_concentration().variable); - read_eigen_matrix(group, SOLID_CONC_DSET, - vars->get_solid_concentration(0).variable); - if (vars->component_has_gas(0)) { - read_eigen_matrix(group, PRESSURE_DSET, - vars->get_pressure_main_variables(0).variable); + auto group = main_vars_grp.open_group(raw_database->get_label_component(0)); + uvars(liquid_saturation) = group.read_vector_dataset(LIQUID_SATURATION_DSET); + uvars(water_aqueous_concentration) = group.read_vector_dataset(AQUEOUS_CONC_DSET); + uvarsi(solid_concentration, 0) = group.read_vector_dataset(SOLID_CONC_DSET); + uvars(capillary_pressure) = group.read_vector_dataset(CAP_PRESSURE_DSET); + if (vars->component_has_gas(0)) + { + uvarsi(pressure_main_variables, 0) + = group.read_vector_dataset(PRESSURE_DSET); } - read_eigen_matrix(group, CAP_PRESSURE_DSET, - vars->get_capillary_pressure().variable); } for (auto component: raw_database->range_aqueous_component()) { - auto group = main_vars_grp.openGroup( + auto group = main_vars_grp.open_group( raw_database->get_label_component(component)); - read_eigen_matrix(group, AQUEOUS_CONC_DSET, - vars->get_aqueous_concentration(component).variable); - read_eigen_matrix(group, SOLID_CONC_DSET, - vars->get_solid_concentration(component).variable); - if (vars->component_has_gas(component)) { - read_eigen_matrix(group, PRESSURE_DSET, - vars->get_pressure_main_variables(component).variable); + uvarsi(aqueous_concentration, component) = group.read_vector_dataset(AQUEOUS_CONC_DSET); + uvarsi(solid_concentration, component) = group.read_vector_dataset(SOLID_CONC_DSET); + if (vars->component_has_gas(component)) + { + uvarsi(pressure_main_variables, component) + = group.read_vector_dataset(PRESSURE_DSET); } } -} + } void UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::read_chemical_solutions( - const H5::Group& timestep_group, + const hdf5::GroupPath& timestep_group, reactmicp::systems::unsaturated::UnsaturatedVariables* vars ) { - H5::Group chem_group = timestep_group.openGroup(CHEMISTRY_SOLUTION_GRPNAME); + hdf5::Group chem_group = timestep_group.open_group(CHEMISTRY_SOLUTION_GRPNAME); database::RawDatabasePtr raw_database = vars->get_database(); mesh::Mesh1DPtr the_mesh = vars->get_mesh(); for (auto node: the_mesh->range_nodes()) { - auto group = chem_group.openGroup(std::to_string(node)); - read_adimensional_system_solution(group, vars->get_adim_solution(node)); + vars->set_adim_solution(node, + read_adimensional_system_solution(chem_group, std::to_string(node))); } } // time to get ugly and use macro // this is simply a quick wrapper that call the read_eigen_matrix function // It avoids all the boilerplate #define read_transport_property(name, var) \ - read_eigen_matrix(trans_group, name, vars->get_##var ().variable) + uvars(var) = trans_group.read_vector_dataset(name); void UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::read_transport_properties( - const H5::Group& timestep_group, + const hdf5::GroupPath& timestep_group, reactmicp::systems::unsaturated::UnsaturatedVariables* vars ) { - H5::Group trans_group = timestep_group.openGroup(TRANSPORT_PROPERTIES_GRPNAME); + hdf5::Group trans_group = timestep_group.open_group(TRANSPORT_PROPERTIES_GRPNAME); read_transport_property(POROSITY_DSET, porosity); read_transport_property(LIQUID_DIFFUSIVITY_DSET, liquid_diffusivity); read_transport_property(REL_LIQUID_DIFFUSIVITY_DSET, relative_liquid_diffusivity); read_transport_property(LIQUID_PERMEABILITY_DSET, liquid_permeability); read_transport_property(REL_LIQUID_PERMEABILITY_DSET, relative_liquid_permeability); read_transport_property(RESISTANCE_GAS_DIFFUSIVITY_DSET, resistance_gas_diffusivity); read_transport_property(REL_GAS_DIFFUSIVITY_DSET, relative_gas_diffusivity); } #undef read_transport_property // we remove the ugly macro, no one saws anything :) +#undef uvars +#undef uvarsi + + +hdf5::Group +UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::open_timestep( + const scalar_t& timestep + ) +{ + auto str_timestep = m_timesteps.get_string(timestep); + return m_file.open_group(str_timestep); +} +hdf5::Group +UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::open_timestep( + const std::string& timestep + ) +{ + return m_file.open_group(timestep); +} + +hdf5::Group +UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::open_main_variables( + const std::string& timestep + ) +{ + auto timestep_grp = open_timestep(timestep); + return timestep_grp.open_group(MAIN_VARIABLES_GRPNAME); +} + +hdf5::Group +UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::open_main_variables_component( + const std::string& timestep, + const std::string& component + ) +{ + auto main_vars_grp = open_main_variables(timestep); + return main_vars_grp.open_group(component); +} + + +//! \brief Return a main variable against x +Vector +UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::main_variable_vs_x( + const std::string& timestep, + const std::string& component, + const std::string& variable + ) +{ + auto comp_grp = open_main_variables_component(timestep, component); + return comp_grp.read_vector_dataset(variable); +} + +//! \brief Return a main variable against t +Vector +UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::main_variable_vs_t( + index_t node, + const std::string& component, + const std::string& variable + ) +{ + Vector data(m_timesteps.size()); + index_t cnt = 0; + for (auto it: m_timesteps) + { + Vector timestep_vector = main_variable_vs_x(it.second, component, variable); + data(cnt) = timestep_vector(node); + ++cnt; + } + return data; +} +hdf5::Group +UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::open_transport_variables( + const std::string& timestep + ) +{ + auto timestep_grp = open_timestep(timestep); + return timestep_grp.open_group(TRANSPORT_PROPERTIES_GRPNAME); +} + +Vector +UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::transport_variable_vs_x( + const std::string& timestep, + const std::string& variable + ) +{ + auto grp = open_transport_variables(timestep); + return grp.read_vector_dataset(variable); +} + +Vector +UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::transport_variable_vs_t( + index_t node, + const std::string& variable + ) +{ + Vector data(m_timesteps.size()); + index_t cnt = 0; + for (auto it: m_timesteps) + { + Vector timestep_vector = transport_variable_vs_x(it.second, variable); + data(cnt) = timestep_vector(node); + ++cnt; + } + return data; +} + +hdf5::Group +UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::open_chemistry_solutions( + const std::string& timestep + ) +{ + auto timestep_grp = open_timestep(timestep); + return timestep_grp.open_group(CHEMISTRY_SOLUTION_GRPNAME); +} + + +AdimensionalSystemSolution +UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::get_adim_solution( + scalar_t timestep, + index_t node + ) +{ + auto str_timestep = m_timesteps.get_string(timestep); + return get_adim_solution(str_timestep, node); +} + +AdimensionalSystemSolution +UnsaturatedHDF5Reader::UnsaturatedHDF5ReaderImpl::get_adim_solution( + std::string timestep, + index_t node + ) +{ + auto chem_grp = open_chemistry_solutions(timestep); + return read_adimensional_system_solution(chem_grp, std::to_string(node)); +} + } //end namespace io } //end namespace specmicp diff --git a/src/reactmicp/io/hdf5_unsaturated.hpp b/src/reactmicp/io/hdf5_unsaturated.hpp index b9de30b..23dd47f 100644 --- a/src/reactmicp/io/hdf5_unsaturated.hpp +++ b/src/reactmicp/io/hdf5_unsaturated.hpp @@ -1,123 +1,176 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. 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. * ============================================================================= */ #ifndef SPECMICP_IO_REACTMICP_HDF5_UNSATURATED_HPP #define SPECMICP_IO_REACTMICP_HDF5_UNSATURATED_HPP //! \file reactmicp/io/hdf5_unsaturated.hpp //! \brief HDF5 formatter for the unsaturated system #include "specmicp_common/types.hpp" +#include "specmicp_common/pimpl_ptr.hpp" #include namespace specmicp { +struct AdimensionalSystemSolution; + namespace units { struct UnitsSet; } //end namespace units namespace reactmicp { namespace systems { namespace unsaturated { class UnsaturatedVariables; } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp namespace io { class HDF5File; +class HDF5Timesteps; //! \brief Save the unsaturated set of variables in HDF5 format //! class SPECMICP_DLL_PUBLIC UnsaturatedHDF5Saver { using UnsaturatedVariables = reactmicp::systems::unsaturated::UnsaturatedVariables; public: UnsaturatedHDF5Saver( const std::string& filename, std::shared_ptr vars, const units::UnitsSet& the_units ); ~UnsaturatedHDF5Saver(); //! \brief Save the variables //! //! All the values will be saved in a new group 'section/name' void save_timestep(scalar_t timestep); private: struct SPECMICP_DLL_LOCAL UnsaturatedHDF5SaverImpl; //! \brief implementation details //! \internal std::unique_ptr m_impl; // delete copy and assignement operator, no need for them UnsaturatedHDF5Saver(const UnsaturatedHDF5Saver& other) = delete; UnsaturatedHDF5Saver& operator= (const UnsaturatedHDF5Saver& other) = delete; }; //! \brief Read saved data from an HDF5 file class SPECMICP_DLL_PUBLIC UnsaturatedHDF5Reader { public: //UnsaturatedHDF5Reader(const HDF5File& file); UnsaturatedHDF5Reader(const std::string& filepath); + ~UnsaturatedHDF5Reader(); + + const HDF5Timesteps& get_timesteps() const; + + units::UnitsSet get_units(); + + AdimensionalSystemSolution get_adim_solution( + std::string timestep, + index_t node + ); + AdimensionalSystemSolution get_adim_solution( + scalar_t timestep, + index_t node + ); + + + //! \brief Return a main variable against x + Vector main_variable_vs_x( + scalar_t timestep, + const std::string& component, + const std::string& variable + ); + //! \brief Return a main variable against x + Vector main_variable_vs_x( + const std::string& timestep, + const std::string& component, + const std::string& variable + ); + //! \brief Return a main variable against t + Vector main_variable_vs_t( + index_t node, + const std::string& component, + const std::string& variable + ); + //! \brief Return a transport variable against x + Vector transport_variable_vs_x( + const std::string& timestep, + const std::string& variable + ); + //! \brief Return a transport variable against x + Vector transport_variable_vs_x( + scalar_t timestep, + const std::string& variable + ); + //! \brief Return a transport variable against t + Vector transport_variable_vs_t( + index_t node, + const std::string& variable + ); //! \brief Initialize variables from a saved timestep //! //! The Mesh and the database must already be ok void initialize_variables( scalar_t timestep, reactmicp::systems::unsaturated::UnsaturatedVariables* vars ); private: struct SPECMICP_DLL_LOCAL UnsaturatedHDF5ReaderImpl; - std::unique_ptr m_impl; + utils::pimpl_ptr m_impl; // delete copy and assignement operator UnsaturatedHDF5Reader(const UnsaturatedHDF5Reader& other) = delete; UnsaturatedHDF5Reader& operator= (const UnsaturatedHDF5Reader& other) = delete; }; } //end namespace io } //end namespace specmicp #endif // SPECMICP_IO_REACTMICP_HDF5_UNSATURATED_HPP diff --git a/tests/reactmicp/CMakeLists.txt b/tests/reactmicp/CMakeLists.txt index 345c7dc..164a8f5 100644 --- a/tests/reactmicp/CMakeLists.txt +++ b/tests/reactmicp/CMakeLists.txt @@ -1,65 +1,66 @@ # Reactmicp : Reactive transport solver # ------------------------------------- set(test_reactmicp_files solver/test_reactive_transport_solver.cpp solver/test_coupling.cpp ) add_catch_test(NAME reactmicp SOURCES ${test_reactmicp_files} LINK_LIBRARIES ${REACTMICP_STATIC_LIBS}) # Saturated diffusion using new reactive transport solver # ------------------------------------------------------- set(test_reactmicp_saturated_files systems/saturated_react/test_reactmicp_saturated_react.cpp systems/saturated_react/speciation_system.cpp systems/saturated_react/variables.cpp systems/saturated_react/transport_program.cpp systems/saturated_react/equilibrium_stagger.cpp ) add_catch_test( NAME reactmicp_saturated_react SOURCES ${test_reactmicp_saturated_files} LINK_LIBRARIES ${REACTMICP_STATIC_LIBS} ) # Unsaturated reactive transport system # ------------------------------------- set(UNSATURATED_DIR systems/unsaturated) set(test_reactmicp_unsaturated_files ${UNSATURATED_DIR}/test_reactmicp_unsaturated.cpp ${UNSATURATED_DIR}/aqueous_equation.cpp ${UNSATURATED_DIR}/aqueous_pressure_equation.cpp ${UNSATURATED_DIR}/boundary_conditions.cpp ${UNSATURATED_DIR}/configuration.cpp ${UNSATURATED_DIR}/equilibrium_stagger.cpp + ${UNSATURATED_DIR}/hdf5_unsaturated.cpp ${UNSATURATED_DIR}/pressure_equation.cpp ${UNSATURATED_DIR}/saturation_equation.cpp ${UNSATURATED_DIR}/saturation_pressure_equation.cpp ${UNSATURATED_DIR}/transport_stagger.cpp ${UNSATURATED_DIR}/variables.cpp ) set(TEST_CEMDATA_PATH \"../../data/cemdata.yaml\") set_source_files_properties( ${test_reactmicp_unsaturated_files} PROPERTIES COMPILE_DEFINITIONS "TEST_CEMDATA_PATH=${TEST_CEMDATA_PATH}" ) add_catch_test( NAME reactmicp_unsaturated SOURCES ${test_reactmicp_unsaturated_files} LINK_LIBRARIES ${REACTMICP_STATIC_LIBS} ) # unsaturated binary # ------------------ add_subdirectory( binary/unsaturated ) diff --git a/tests/reactmicp/systems/unsaturated/hdf5_unsaturated.cpp b/tests/reactmicp/systems/unsaturated/hdf5_unsaturated.cpp new file mode 100644 index 0000000..ea9daba --- /dev/null +++ b/tests/reactmicp/systems/unsaturated/hdf5_unsaturated.cpp @@ -0,0 +1,152 @@ +#include "catch.hpp" + +#include "reactmicp/io/hdf5_unsaturated.hpp" + +#include "specmicp/adimensional/adimensional_system_solution.hpp" +#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" + +#include "specmicp_database/database.hpp" + +#include "specmicp/problem_solver/reactant_box.hpp" +#include "specmicp/adimensional/adimensional_system_solver.hpp" + + +#include "reactmicp/systems/unsaturated/variables_interface.hpp" + +#include "dfpm/mesh.hpp" + +#include + +#define PATH_TEST_FILES "test_hdf5_unsaturated" + +using namespace specmicp; +using namespace specmicp::reactmicp::systems::unsaturated; + +static specmicp::database::RawDatabasePtr get_database() +{ + static database::RawDatabasePtr raw_data {nullptr}; + if (raw_data == nullptr) + { + specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); + thedatabase.keep_only_components( + {"H2O", "H[+]", "Ca[2+]", "Si(OH)4", "HCO3[-]"} + ); + std::map swap = {{"HCO3[-]", "CO2"},}; + thedatabase.swap_components(swap); + raw_data = thedatabase.get_database(); + raw_data->freeze_db(); + } + return raw_data; +} + +static AdimensionalSystemSolution +get_init_solution( + database::RawDatabasePtr the_database, + units::UnitsSet the_units + ) +{ + + specmicp::ReactantBox reactant_box(the_database, the_units); + reactant_box.set_solution(500, "L"); + reactant_box.add_solid_phase("C3S", 350, "L"); + reactant_box.add_solid_phase("C2S", 50, "L"); + reactant_box.add_aqueous_species("CO2", 32.72, "mol/kg"); + reactant_box.set_charge_keeper("HO[-]"); + + + auto constraints = reactant_box.get_constraints(true); + + AdimensionalSystemSolver adim_solver(the_database, constraints); + adim_solver.get_options().solver_options.set_tolerance(1e-14); + adim_solver.get_options().units_set = the_units; + + Vector x; + adim_solver.initialise_variables(x, 0.8, -3); + micpsolver::MiCPPerformance perf = adim_solver.solve(x); + specmicp_assert(perf.return_code > micpsolver::MiCPSolverReturnCode::Success); + + return adim_solver.get_raw_solution(x); +} + + + +TEST_CASE("hdf5 output", "[hdf5]") { + + std::string filename = "test_hdf5_unsaturated.hdf5"; + + database::RawDatabasePtr raw_data = get_database(); + + mesh::Uniform1DMeshGeometry mesh_geom; + mesh_geom.nb_nodes = 10; + mesh_geom.dx = 1; + mesh_geom.section = 1; + mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d(mesh_geom); + + + units::UnitsSet the_units; + the_units.length = units::LengthUnit::centimeter; + + specmicp::AdimensionalSystemSolution init_sol + = get_init_solution(raw_data, the_units); + + specmicp::AdimensionalSystemSolutionExtractor extr(init_sol, raw_data, the_units); + + SECTION("Writer") { + + + VariablesInterface vars_init(the_mesh, raw_data, + {raw_data->get_id_component_from_element("C")}, the_units); + + for (auto node: the_mesh->range_nodes()) + { + vars_init.initialize_variables(node, extr); + } + + auto vars = vars_init.get_variables(); + io::UnsaturatedHDF5Saver saver(filename, vars, the_units); + saver.save_timestep(10.0); + } + + + + SECTION("Reader") { + + specmicp::io::UnsaturatedHDF5Reader reader(filename); + + // units + units::UnitsSet read_unit = reader.get_units(); + CHECK(read_unit.length == the_units.length); + CHECK(read_unit.mass == the_units.mass); + CHECK(read_unit.quantity == the_units.quantity); + + // adimensional solution + specmicp::AdimensionalSystemSolution sol = reader.get_adim_solution(10.0, 1); + + REQUIRE(sol.is_valid); + for (auto r: range(sol.main_variables.rows())) + { + if (std::isfinite(sol.main_variables(r))) + { + CHECK(sol.main_variables(r) == Approx(init_sol.main_variables(r)).epsilon(1e-10)); + } + } + CHECK(sol.log_gamma.norm() == Approx(init_sol.log_gamma.norm()).epsilon(1e-10)); + CHECK(sol.gas_fugacities.norm() == Approx(init_sol.gas_fugacities.norm()).epsilon(1e-10)); + CHECK(sol.secondary_molalities.norm() == Approx(init_sol.secondary_molalities.norm()).epsilon(1e-10)); + + + // main variables + Vector liq_sat = reader.main_variable_vs_x(10.0, "H2O", "liquid_saturation"); + CHECK(liq_sat.rows() == the_mesh->nb_nodes()); + CHECK(liq_sat(1) == Approx(extr.saturation_water())); + + Vector ca_s = reader.main_variable_vs_x(10.0, "Ca[2+]", "solid_concentration"); + CHECK(ca_s.rows() == the_mesh->nb_nodes()); + CHECK(ca_s(1) == Approx(extr.total_solid_concentration(raw_data->get_id_component("Ca[2+]")))); + + + + + + } +}