diff --git a/examples/reactmicp/unsaturated/drying.cpp b/examples/reactmicp/unsaturated/drying.cpp index a9abacc..0f1d34b 100644 --- a/examples/reactmicp/unsaturated/drying.cpp +++ b/examples/reactmicp/unsaturated/drying.cpp @@ -1,578 +1,580 @@ /* ============================================================================= 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. * ============================================================================= */ // ====================================================== // // // // Drying example // // ============== // // // // ====================================================== // // // // This is a simple example to show how the unsaturated // system of ReactMiCP is able to solve a complex drying // equation. // // It implements a special chemistry stagger (DryingStagger) // which only computes the vapor pressure. // Macro to define the problem // =========================== #define NB_NODES 50 #define DX 0.05*1e-2 // cm #define CROSS_SECTION 5.0*1e-4 // cm^2 #define GAMMA_WATER_GLASS 0.1 // N/m #define GAMMA_WATER_AIR 71.97e-3 // N/m = kg/s^2 #define R_SMALL_BEAD 0.015*1e-2 // cm #define R_BIG_BEAD 0.04*1e-2 // cm #define M_v 18.015e-3 // kg/mol // Includes // ======== #include "specmicp_database/database.hpp" #include "reactmicp/systems/unsaturated/variables.hpp" #include "reactmicp/systems/unsaturated/variables_sub.hpp" #include "reactmicp/systems/unsaturated/transport_stagger.hpp" #include "reactmicp/solver/staggers_base/chemistry_stagger_base.hpp" #include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" +#include "reactmicp/solver/reactive_transport_solver.hpp" +#include "reactmicp/solver/reactive_transport_solver_structs.hpp" #include "reactmicp/solver/runner.hpp" #include "dfpm/meshes/generic_mesh1d.hpp" #include "specmicp_common/physics/laws.hpp" #include "specmicp_common/physics/units.hpp" #include "specmicp_common/log.hpp" #include "specmicp_common/io/csv_formatter.hpp" #include #include // Namespace // ========= using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems; using namespace specmicp::reactmicp::systems::unsaturated; // Declarations // ============ // Constants // ========== static const scalar_t rho_l_si = constants::water_density_25; static const scalar_t rho_l = 1e-6*constants::water_density_25; static const scalar_t cw_tilde_si = rho_l_si / M_v; static const scalar_t cw_tilde = 1e-6*cw_tilde_si; int main(); database::RawDatabasePtr get_database(); mesh::Mesh1DPtr get_mesh(); scalar_t leverett_function(scalar_t s_eff); scalar_t big_bead_cap_pressure(scalar_t saturation); scalar_t small_bead_cap_pressure(scalar_t saturation); scalar_t vapor_pressure(scalar_t pc); // Drying Stagger // ============== // // This stagger computes the exchange between water vapor and liquid water class DryingStagger: public solver::ChemistryStaggerBase { public: DryingStagger() {} virtual void initialize(solver::VariablesBase * const var) override {} //! \brief Initialize the stagger at the beginning of an iteration //! //! This is where the predictor can be saved, the first trivial iteration done, ... //! //! \param dt the duration of the timestep //! \param var a shared_ptr to the variables virtual void initialize_timestep( scalar_t dt, solver::VariablesBase * const var ) override { m_dt = dt; } //! \brief Solve the equation for the timestep //! //! \param var a shared_ptr to the variables virtual solver::StaggerReturnCode restart_timestep( solver::VariablesBase * const var) override; void compute_one_node(index_t node, UnsaturatedVariables * const vars); private: scalar_t m_dt {-1.0}; }; // Upscaling Stagger // ================= // // Define how we compute capillary pressure, vapor pressure, ... class UpscalingDryingStagger: public solver::UpscalingStaggerBase { public: UpscalingDryingStagger(std::vector& is_big_bead): m_is_big_bead(is_big_bead) {} //! \brief Initialize the stagger at the beginning of the computation //! //! \param var a shared_ptr to the variables virtual void initialize(solver::VariablesBase * const var) override; //! \brief Initialize the stagger at the beginning of an iteration //! //! This is where the predictor can be saved, the first trivial iteration done, ... //! //! \param dt the duration of the timestep //! \param var a shared_ptr to the variables virtual void initialize_timestep( scalar_t dt, solver::VariablesBase * const var ) override {} //! \brief Solve the equation for the timestep //! //! \param var a shared_ptr to the variables virtual solver::StaggerReturnCode restart_timestep( solver::VariablesBase * const var) override { return solver::StaggerReturnCode::ResidualMinimized; } scalar_t capillary_pressure_model(index_t node, scalar_t saturation); scalar_t vapor_pressure_model(index_t node, scalar_t saturation); scalar_t relative_liquid_permeability_model(index_t node, scalar_t saturation); scalar_t relative_liquid_diffusivity_model(index_t node, scalar_t saturation); scalar_t relative_gas_diffusivity_model(index_t node, scalar_t saturation); private: std::vector m_is_big_bead; }; // Implementation // ============== scalar_t UpscalingDryingStagger::capillary_pressure_model(index_t node, scalar_t saturation) { scalar_t pc = 0.0; if (m_is_big_bead[node]) { pc = big_bead_cap_pressure(saturation); } else pc = small_bead_cap_pressure(saturation); return pc; } scalar_t UpscalingDryingStagger::vapor_pressure_model(index_t node, scalar_t saturation) { return vapor_pressure(capillary_pressure_model(node, saturation)); } scalar_t UpscalingDryingStagger::relative_liquid_permeability_model(index_t node, scalar_t saturation) { return std::pow(saturation, 3); } scalar_t UpscalingDryingStagger::relative_liquid_diffusivity_model(index_t node, scalar_t saturation) { return saturation; } scalar_t UpscalingDryingStagger::relative_gas_diffusivity_model(index_t node, scalar_t saturation) { return (1.0 - saturation); } void UpscalingDryingStagger::initialize(solver::VariablesBase * const var) { UnsaturatedVariables * const true_vars = static_cast(var); true_vars->set_capillary_pressure_model( [this](index_t node, scalar_t sat) -> scalar_t { return this->capillary_pressure_model(node, sat);}); true_vars->set_vapor_pressure_model(std::bind(std::mem_fn( &UpscalingDryingStagger::vapor_pressure_model ), this, std::placeholders::_1, std::placeholders::_2)); true_vars->set_relative_liquid_permeability_model(std::bind(std::mem_fn( &UpscalingDryingStagger::relative_liquid_permeability_model ), this, std::placeholders::_1, std::placeholders::_2)); true_vars->set_relative_liquid_diffusivity_model(std::bind(std::mem_fn( &UpscalingDryingStagger::relative_liquid_diffusivity_model ), this, std::placeholders::_1, std::placeholders::_2)); true_vars->set_relative_gas_diffusivity_model(std::bind(std::mem_fn( &UpscalingDryingStagger::relative_gas_diffusivity_model ), this, std::placeholders::_1, std::placeholders::_2)); MainVariable& vapor_pressure = true_vars->get_pressure_main_variables(0); MainVariable& saturation = true_vars->get_liquid_saturation(); SecondaryTransientVariable& porosity = true_vars->get_porosity(); SecondaryVariable& liquid_diffusivity = true_vars->get_liquid_diffusivity(); SecondaryVariable& liquid_permeability = true_vars->get_liquid_permeability(); SecondaryVariable& resistance_gas_diffusivity = true_vars->get_resistance_gas_diffusivity(); true_vars->get_binary_gas_diffusivity(0) = 0.282*1e-4; porosity(0) = 1; resistance_gas_diffusivity(0) = 1.0; for (index_t node=1; nodeget_mesh()->nb_nodes(); ++node) { vapor_pressure(node) = vapor_pressure_model(node, saturation(node)); if (m_is_big_bead[node]) { porosity(node) = 0.371; liquid_diffusivity(node) = 0; liquid_permeability(node) = 3.52e-11; resistance_gas_diffusivity(node) = 2*0.371/(3-0.371); } else { porosity(node) = 0.387; liquid_diffusivity(node) = 0; liquid_permeability(node) = 8.41e-12; resistance_gas_diffusivity(node) = 2*0.387/(3-0.387); } } porosity.predictor = porosity.variable; true_vars->set_relative_variables(); } solver::StaggerReturnCode DryingStagger::restart_timestep( solver::VariablesBase * const var) { UnsaturatedVariables * const true_vars = static_cast(var); for (index_t node=1; nodeget_mesh()->nb_nodes(); ++node) { compute_one_node(node, true_vars); } return solver::StaggerReturnCode::ResidualMinimized; } void DryingStagger::compute_one_node( index_t node, UnsaturatedVariables * const vars) { MainVariable& satvars = vars->get_liquid_saturation(); MainVariable& presvars = vars->get_pressure_main_variables(0); scalar_t rt = vars->get_rt(); scalar_t saturation = satvars(node); scalar_t pressure = presvars(node); scalar_t porosity = vars->get_porosity()(node); auto vapor_pressure_f = vars->get_vapor_pressure_model(); const scalar_t tot_conc = cw_tilde_si*saturation + (1.0-saturation)*pressure/rt; pressure = vapor_pressure_f(node, saturation); scalar_t cw_hat = pressure/rt; saturation = (tot_conc - cw_hat) / (cw_tilde_si-cw_hat); pressure = vapor_pressure_f(node, saturation); cw_hat = pressure/rt; scalar_t res = tot_conc - cw_tilde_si*saturation - (1.0-saturation)*cw_hat; while (std::abs(res)/tot_conc > 1e-12) { saturation = (tot_conc - cw_hat) / (cw_tilde_si-cw_hat); pressure = vapor_pressure_f(node, saturation); cw_hat = pressure/rt; res = tot_conc - cw_tilde_si*saturation - (1.0-saturation)*cw_hat; } satvars(node) = saturation; satvars.velocity(node) = (saturation - satvars.predictor(node))/m_dt; presvars(node) = pressure; presvars.velocity(node) = (pressure - presvars.predictor(node))/m_dt; presvars.chemistry_rate(node) = - porosity/(rt*m_dt) * ( (pressure*(1.0-saturation)) - (presvars.predictor(node)*(1.0-satvars.predictor(node))) ) + presvars.transport_fluxes(node); } database::RawDatabasePtr get_database() { specmicp::database::Database thedatabase(CEMDATA_PATH); std::vector to_keep = {"H2O", "H[+]"}; thedatabase.keep_only_components(to_keep); thedatabase.remove_gas_phases(); database::RawDatabasePtr raw_data = thedatabase.get_database(); raw_data->freeze_db(); return raw_data; } mesh::Mesh1DPtr get_mesh() { mesh::Uniform1DMeshGeometry geom; geom.dx = DX; geom.nb_nodes = NB_NODES; geom.section = CROSS_SECTION; return mesh::uniform_mesh1d(geom); } scalar_t big_bead_cap_pressure(scalar_t saturation) { return GAMMA_WATER_AIR/std::sqrt(3.52e-11/0.371)*leverett_function(saturation); } scalar_t small_bead_cap_pressure(scalar_t saturation) { return GAMMA_WATER_AIR/std::sqrt(8.41e-12/0.387)*leverett_function(saturation); } scalar_t vapor_pressure(scalar_t pc) { return 3e3*std::exp(-M_v * pc / (rho_l_si * (8.314*(25+273.15)))); } scalar_t leverett_function(scalar_t s_eff) { return 0.325*std::pow((1.0/s_eff - 1), 0.217); } class Printer { public: Printer(UnsaturatedVariables* vars, std::size_t reserve_size); void register_timestep(scalar_t time, solver::VariablesBasePtr base_vars); void print(const std::string& name); private: mesh::Mesh1DPtr m_mesh; std::vector> m_pressure; std::vector> m_saturation; std::vector m_timestep; }; Printer::Printer(UnsaturatedVariables *vars, std::size_t reserve_size): m_mesh(vars->get_mesh()), m_pressure(m_mesh->nb_nodes()), m_saturation(m_mesh->nb_nodes()) { m_timestep.reserve(reserve_size); for (auto node: m_mesh->range_nodes()) { m_pressure[node].reserve(reserve_size); m_saturation[node].reserve(reserve_size); } } void Printer::register_timestep(scalar_t time, solver::VariablesBasePtr base_vars) { UnsaturatedVariables* vars = static_cast(base_vars.get()); m_timestep.push_back(time); MainVariable& saturation = vars->get_liquid_saturation(); MainVariable& pressure = vars->get_pressure_main_variables(0); for (auto node: m_mesh->range_nodes()) { m_pressure[node].push_back(pressure(node)); m_saturation[node].push_back(saturation(node)); } } void Printer::print(const std::string& name) { io::CSVFile pressure_file(name+"_pressure.dat"); io::CSVFile saturation_file(name+"_saturation.dat"); pressure_file << "Node"; saturation_file << "Node"; for (auto& time: m_timestep) { pressure_file.separator(); pressure_file << time; saturation_file.separator(); saturation_file << time; } pressure_file.eol(); saturation_file.eol(); for (auto node: m_mesh->range_nodes()) { pressure_file << m_mesh->get_position(node); saturation_file << m_mesh->get_position(node); for (std::size_t ind=0; ind has_gas = {true, false, false}; std::vector is_big_bead(the_mesh->nb_nodes(), true); //for (int node=25; node( the_mesh, raw_data, has_gas, units::SI_units); auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); bcs->add_gas_node(0); bcs->fork_constraint(0, "gas_node"); // variables initialisation MainVariable& saturation = vars->get_liquid_saturation(); MainVariable& vapor_pressure = vars->get_pressure_main_variables(0); SecondaryTransientVariable& water_aqueous_concentration = vars->get_water_aqueous_concentration(); water_aqueous_concentration.predictor = water_aqueous_concentration.variable; SecondaryVariable& rel_gas_d = vars->get_relative_gas_diffusivity(); rel_gas_d(0) = 1.0; vapor_pressure(0) = 100; for (index_t node = 1; node < the_mesh->nb_nodes(); ++node) { saturation(node) = 0.9; } water_aqueous_concentration.set_constant(cw_tilde_si); std::shared_ptr upscaling_stagger = std::make_shared(is_big_bead); upscaling_stagger->initialize(vars.get()); std::shared_ptr drying_stagger = std::make_shared(); // init vapor pressure drying_stagger->initialize_timestep(1.0, vars.get()); drying_stagger->restart_timestep(vars.get()); std::shared_ptr transport_stagger = std::make_shared(vars, bcs); solver::ReactiveTransportSolver react_solver(transport_stagger, drying_stagger, upscaling_stagger); solver::ReactiveTransportOptions& opts = react_solver.get_options(); opts.residuals_tolerance = 1e-2; opts.good_enough_tolerance = 0.99; vapor_pressure.chemistry_rate.setZero(); //for (auto it=0; it<100; ++it) //{ // auto retcode = react_solver.solve_timestep(0.1, vars); // std::cout << (int) retcode << std::endl; //} //std::cout << saturation.variable << std::endl; //std::cout << vapor_pressure.variable << std::endl; //std::cout << vars->get_gas_diffusivity().variable << std::endl; //std::cout << vars->get_relative_gas_diffusivity().variable << std::endl; int nb_hours = 9; Printer printer(vars.get(), 4*nb_hours); solver::SimulationInformation simul_info("drying", 900); solver::ReactiveTransportRunner runner(react_solver, 0.01, 10, simul_info); solver::output_f out_func = std::bind(&Printer::register_timestep, &printer, std::placeholders::_1, std::placeholders::_2); runner.set_output_policy(out_func); runner.run_until(nb_hours*3600, vars); printer.print("out_drying"); std::cout << saturation.variable << std::endl; std::cout << vapor_pressure.variable << std::endl; } diff --git a/src/bin/reactmicp/unsaturated.cpp b/src/bin/reactmicp/unsaturated.cpp index 211bac8..361dc85 100644 --- a/src/bin/reactmicp/unsaturated.cpp +++ b/src/bin/reactmicp/unsaturated.cpp @@ -1,856 +1,854 @@ /* ============================================================================= 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. * ============================================================================= */ //! \file bin/reactmicp/unsaturated.cpp //! \brief Unsaturated #include "specmicp_common/types.hpp" #include "specmicp_common/cli/parser.hpp" #include "specmicp_common/plugins/plugin_manager.hpp" #include "specmicp_common/io/yaml.hpp" #include "specmicp_common/io/safe_config.hpp" #include "specmicp_common/compat.hpp" #include "specmicp_common/log.hpp" #include "specmicp_common/filesystem.hpp" #include "specmicp_common/string_algorithms.hpp" #include "specmicp_common/openmp_support.hpp" #include "specmicp_common/physics/io/configuration.hpp" #include "specmicp_database/io/configuration.hpp" #include "dfpm/io/configuration.hpp" #include "specmicp_common/io/configuration.hpp" #include "specmicp/io/print.hpp" #include "reactmicp/reactmicp_unsaturated.hpp" #include "reactmicp/systems/unsaturated/variables.hpp" #include "dfpm/solver/parabolic_structs.hpp" #include "reactmicp/io/hdf5_unsaturated.hpp" #include "specmicp_common/config.h" #include #include #include #include #include "specmicp_common/io/config_yaml_sections.h" #include "specmicp_common/plugins/plugin_modules_names.h" using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems::unsaturated; const char* welcome = R"(=============================== = Reactmicp - unsaturated = =============================== (C) copyright 2014-2016 F.Georget, Princeton University version: not good enough yet )"; int main (int argc, char* argv[]); plugins::PluginManager& initialize_plugin_manager( io::YAMLConfigHandle& config ); SimulationData initialize_simul_data( io::YAMLConfigHandle& safe_config, const std::vector& db_dirs ); std::shared_ptr initialize_variables( std::shared_ptr raw_data, std::shared_ptr the_mesh, const units::UnitsSet& the_units, io::YAMLConfigHandle&& configuration ); struct Staggers { std::shared_ptr upscaling; std::shared_ptr chemistry; std::shared_ptr transport; }; Staggers initialize_staggers( SimulationData& simul_data, io::YAMLConfigHandle&& configuration ); std::shared_ptr initialize_upscaling_stagger( SimulationData& simul_data, io::YAMLConfigHandle& conf_staggers ); std::shared_ptr initialize_chemistry_stagger( SimulationData& simul_data, io::YAMLConfigHandle& conf_staggers ); std::shared_ptr initialize_transport_stagger( SimulationData& simul_data, io::YAMLConfigHandle& conf_staggers ); std::string initialize_output( SimulationData& simul_data, const std::string& work_dir, io::YAMLConfigHandle&& conf ); void print_db_to_conf_log(database::DataContainer* raw_data); int main(int argc, char *argv[]) { std::cout << welcome << std::endl; init_logger(&std::cout, logger::Warning); // CLI options cli::CommandLineParser parser; parser.add_option('i', "input", cli::ValueType::string, "input_file (yaml format)"); parser.add_option('d', "database_dir", std::string(""), "directory where database are stored"); parser.add_option('w', "workdir", utils::get_current_directory(), "working directory"); parser.set_help_message("ReactMiCP Unsaturated System"); const auto ret_parse = parser.parse(argc, argv); if (ret_parse > 0) { return EXIT_SUCCESS; } - std::string input_filepath = parser.get_option("input"); std::string database_dir = parser.get_option("database_dir"); std::string working_dir = parser.get_option("workdir"); std::vector db_dirs; const auto env = utils::get_env("SPECMICP_DATABASE_PATH"); if (not env.empty()) { db_dirs = utils::split(env, ';'); } if (not database_dir.empty()) {db_dirs.push_back(database_dir);} if (not working_dir.empty()) {db_dirs.push_back(working_dir);} // Configuration File // ------------------ auto config = io::YAMLConfigFile::make(input_filepath); std::unique_ptr out_conf_logger; if (config->has_section(SPC_CF_S_CONF_LOGS)) { out_conf_logger = io::configure_conf_log( config->get_section(SPC_CF_S_CONF_LOGS), working_dir ); } SPC_CONF_LOG << "Running ReactMiCP unsaturated driver"; std::unique_ptr out_logger; if (config->has_section(SPC_CF_S_LOGS)) { out_logger = io::configure_log( config->get_section(SPC_CF_S_LOGS), working_dir ); } initialize_plugin_manager(*(config.get())); // openmp initialization // --------------------- utils::initialize_parallel(); // Simulation Data // --------------- SimulationData simul_data = initialize_simul_data(*config.get(), db_dirs); // Staggers // ------- Staggers staggers = initialize_staggers( simul_data, config->get_section(SPC_CF_S_STAGGERS) ); // Check the variables // ------------------- auto res = VariablesInterface(simul_data.variables).check_variables(); if (res != VariablesValidity::good) { const auto msg = "Variables were not initialized correctly." "Please see log for errors."; ERROR << msg; throw std::invalid_argument(msg); } // Solver // ------ solver::ReactiveTransportSolver react_solver( staggers.transport, staggers.chemistry, staggers.upscaling); io::configure_reactmicp_options( react_solver.get_options(), config->get_section(SPC_CF_S_REACTMICP) ); SPC_CONF_LOG << "\n\n++ ReactMiCP ++\n" << SPC_CONF_LOG_SECTION << "\n ReactMiCP options\n" << SPC_CONF_LOG_HLINE; std::ostringstream msg; io::print_reactmicp_options(msg, react_solver.get_options()); SPC_CONF_LOG << msg.str(); // Initialize output std::unique_ptr saver; solver::SimulationInformation simul_info = io::configure_simulation_information( config->get_section(SPC_CF_S_SIMULINFO)); simul_info.working_dir = working_dir; solver::ReactiveTransportRunner runner(react_solver, 1.0, 10.0, // dummy value, set later simul_info); io::configure_reactmicp_timestepper( runner.get_timestepper_options(), config->get_section(SPC_CF_S_TIMESTEPPER) ); // if (config->has_section(SPC_CF_S_REACTOUTPUT)) { auto filepath = initialize_output( simul_data, working_dir, config->get_section(SPC_CF_S_REACTOUTPUT)); saver = make_unique( filepath, simul_data.variables, simul_data.units); auto* saver_ptr = saver.get(); auto out_pol = [saver_ptr]( scalar_t timestep, solver::VariablesBasePtr _) { saver_ptr->save_timestep(timestep); }; saver->save_timestep(0.0); runner.set_output_policy(out_pol); } // run info auto run_section = config->get_section(SPC_CF_S_RUN); auto run_until = run_section.get_required_attribute( SPC_CF_S_RUN_A_RUNUTNIL); // we destroy config to get information config.reset(nullptr); runner.run_until(run_until, simul_data.variables); return EXIT_SUCCESS; } // ==================== // // Simulation data // // ==================== // SimulationData initialize_simul_data( io::YAMLConfigHandle& safe_config, const std::vector& db_dirs ) { // Package the information SimulationData simul_data; // The units // --------- if (safe_config.has_section(SPC_CF_S_UNITS)) { simul_data.units = io::configure_units( safe_config.get_section(SPC_CF_S_UNITS)); } SPC_CONF_LOG << SPC_CONF_LOG_HLINE << "\nUnits :" << "\n - length : " << io::to_string(simul_data.units.length) << "\n - mass : " << io::to_string(simul_data.units.mass) << "\n - quantity : " << io::to_string(simul_data.units.quantity) << "\n" << SPC_CONF_LOG_HLINE; // The database // ------------ SPC_CONF_LOG << "Database initialization\n"; simul_data.raw_data = io::configure_database( safe_config.get_section(SPC_CF_S_DATABASE), db_dirs ); SPC_CONF_LOG << SPC_CONF_LOG_HLINE; if (not simul_data.raw_data->is_valid()) { CRITICAL << "Invalid database from the configuration."; throw std::runtime_error("Database parsed from the configuration has" " been detected to be invalid. Abort."); } // The mesh // -------- SPC_CONF_LOG << "Mesh initialization\n" << SPC_CONF_LOG_SECTION; simul_data.mesh1d = io::configure_mesh(safe_config.get_section(SPC_CF_S_MESH)); SPC_CONF_LOG << SPC_CONF_LOG_HLINE; // The variables // ------------- SPC_CONF_LOG << "Variables initialization\n" << SPC_CONF_LOG_SECTION; simul_data.variables = initialize_variables( simul_data.raw_data, simul_data.mesh1d, simul_data.units, safe_config.get_section(SPC_CF_S_REACTMICPVARIABLES) ); if (not simul_data.raw_data->is_valid()) { CRITICAL << "Invalid database after variable initialization."; throw std::runtime_error("Database has been detected to be invalid" " after the variables initialization. Abort."); } SPC_CONF_LOG << SPC_CONF_LOG_HLINE; // The database is supposed set after this point ! simul_data.raw_data->freeze_db(); print_db_to_conf_log(simul_data.raw_data.get()); // The boundary conditions // ----------------------- SPC_CONF_LOG << "Boundary conditions initialization\n" << SPC_CONF_LOG_SECTION; simul_data.bcs = io::configure_unsaturated_boundary_conditions( simul_data.mesh1d->nb_nodes(), simul_data.raw_data.get(), safe_config.get_section(SPC_CF_S_BOUNDARYCONDITIONS) ); SPC_CONF_LOG << SPC_CONF_LOG_HLINE; // Simulation data ready return simul_data; } - std::shared_ptr initialize_variables( std::shared_ptr raw_data, std::shared_ptr the_mesh, const units::UnitsSet& the_units, io::YAMLConfigHandle&& configuration ) { auto type = configuration.get_optional_attribute( SPC_CF_S_REACTMICPVARIABLES_A_TYPE, SPC_CF_V_PLUGIN); if (type == SPC_CF_V_DEFAULT) type = SPC_CF_V_PLUGIN; if (type != SPC_CF_V_PLUGIN) { configuration.report_error( io::YAMLConfigError::InvalidArgument, "Only accepted value for attribute '" SPC_CF_S_REACTMICPVARIABLES_A_TYPE "' is '" SPC_CF_V_PLUGIN "'." ); } auto plugin_name = configuration.get_required_attribute( SPC_CF_A_PLUGIN_NAME ); SPC_CONF_LOG << "Plugin name : " << plugin_name; plugins::PluginManager& plugin_manager = plugins::get_plugin_manager(); auto extra_plugin_file = configuration.get_optional_attribute( SPC_CF_A_PLUGIN_FILE, ""); if (extra_plugin_file != "") { SPC_CONF_LOG << "Load plugin file : " << extra_plugin_file; auto retcode = plugin_manager.load_plugin(extra_plugin_file); if (not retcode) { CRITICAL << "Failed to load plugin file : " << extra_plugin_file << std::endl; std::runtime_error("Fail to load plugin file : " + extra_plugin_file + "."); } } auto initializer = plugin_manager.get_object( SPC_PLUGIN_REACTMICP_UNSATURATED_VARIABLES, plugin_name); if (initializer == nullptr) { CRITICAL << "No variable initializer found"; throw std::runtime_error("No variable initializer found"); } return initializer->initialize_variables( raw_data, the_mesh, the_units, configuration); } // ================ // // Plugin Manager // // ================ // plugins::PluginManager& initialize_plugin_manager( io::YAMLConfigHandle& config ) { SPC_CONF_LOG << SPC_CONF_LOG_SECTION << "\nInitialize plugin manager\n" << SPC_CONF_LOG_HLINE; plugins::PluginManager& plugin_manager = plugins::get_plugin_manager(); SPC_CONF_LOG << "Registering module " SPC_PLUGIN_REACTMICP_UNSATURATED_VARIABLES; plugin_manager.register_module(SPC_PLUGIN_REACTMICP_UNSATURATED_VARIABLES, make_unique()); SPC_CONF_LOG << "Registering module " SPC_PLUGIN_REACTMICP_UNSATURATED_UPSCALING; plugin_manager.register_module(SPC_PLUGIN_REACTMICP_UNSATURATED_UPSCALING, make_unique()); if (config.has_section(SPC_CF_S_PLUGINS)) { io::configure_plugin_manager(config.get_section(SPC_CF_S_PLUGINS)); } auto dirs = plugin_manager.get_plugin_directories(); if (not dirs.empty()) { SPC_CONF_LOG << "plugin path : "; for (auto dir: dirs) { SPC_CONF_LOG << " - " << dir; } } SPC_CONF_LOG << SPC_CONF_LOG_HLINE; return plugin_manager; } // ================ // // Staggers // // ================ // Staggers initialize_staggers( SimulationData& simul_data, io::YAMLConfigHandle&& conf_staggers ) { SPC_CONF_LOG << "\n\n" << SPC_CONF_LOG_SECTION << "\n= Staggers initialization =\n" << SPC_CONF_LOG_SECTION; Staggers staggers; auto* var_ptr = simul_data.variables.get(); staggers.upscaling = initialize_upscaling_stagger(simul_data, conf_staggers); staggers.upscaling->initialize(var_ptr); staggers.chemistry = initialize_chemistry_stagger(simul_data, conf_staggers); staggers.chemistry->initialize(var_ptr); staggers.transport = initialize_transport_stagger(simul_data, conf_staggers); staggers.transport->initialize(var_ptr); staggers.upscaling->register_transport_stagger(staggers.transport); staggers.upscaling->register_chemistry_stagger(staggers.chemistry); // check the database if (not simul_data.raw_data->is_valid()) { CRITICAL << "Change in database during stagger initialization !"; throw std::logic_error("Change in the database has been detected" " during the staggers initialization. Abort."); } return staggers; } std::shared_ptr initialize_upscaling_stagger( SimulationData& simul_data, io::YAMLConfigHandle& config_staggers ) { SPC_CONF_LOG << "\n\n --> Upscaling stagger <-- \n" << SPC_CONF_LOG_SECTION; auto conf = config_staggers.get_section(SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER); auto type = conf.get_optional_attribute( SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER_A_TYPE, SPC_CF_V_PLUGIN); // only plugin accepted for now if (type == SPC_CF_V_DEFAULT) type = SPC_CF_V_PLUGIN; if (type != SPC_CF_V_PLUGIN) { conf.report_error( io::YAMLConfigError::InvalidArgument, "Invalid argument for attribute '" SPC_CF_A_TYPE "' only value accepted is '" SPC_CF_V_PLUGIN "'."); } auto plugin_name = conf.get_required_attribute( SPC_CF_A_PLUGIN_NAME); SPC_CONF_LOG << "Use plugin " << plugin_name; plugins::PluginManager& plugin_manager = plugins::get_plugin_manager(); auto extra_plugin_file = conf.get_optional_attribute( SPC_CF_A_PLUGIN_FILE, ""); if (extra_plugin_file != "") { SPC_CONF_LOG << "Load plugin file " << extra_plugin_file; auto retcode = plugin_manager.load_plugin(extra_plugin_file); if (not retcode) { CRITICAL << "Failed to load plugin file : " << extra_plugin_file << std::endl; std::runtime_error("Fail to load plugin file : " + extra_plugin_file + "."); } } auto upscaling_factory = plugin_manager.get_object( SPC_PLUGIN_REACTMICP_UNSATURATED_UPSCALING, plugin_name); if (upscaling_factory == nullptr) { CRITICAL << "No upscaling stagger factory found"; throw std::runtime_error("No upscaling stagger factory found"); } SPC_CONF_LOG << SPC_CONF_LOG_HLINE; return upscaling_factory->get_upscaling_stagger( simul_data, std::move(conf)); } std::shared_ptr initialize_chemistry_stagger( SimulationData& simul_data, io::YAMLConfigHandle& conf_staggers ) { SPC_CONF_LOG << "\n\n --> Chemistry stagger <-- \n" << SPC_CONF_LOG_SECTION; // no conf // ------- if (not conf_staggers.has_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER)) { auto opts = EquilibriumOptionsVector::make( simul_data.mesh1d->nb_nodes()); return EquilibriumStagger::make( simul_data.variables, simul_data.bcs, opts ); } // conf is given // -------------- auto conf_chem = conf_staggers.get_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER); auto type = conf_chem.get_optional_attribute( SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE, SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM); if (type == SPC_CF_V_DEFAULT) type = SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM; if (type !=SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM) { CRITICAL << "Unknown chemistry stagger"; throw std::invalid_argument("Only equilibrium stagger supported."); } SPC_CONF_LOG << "Use equilibrium stagger"; std::shared_ptr opts {nullptr}; if (conf_chem.has_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_SS_EQUILIBRIUM_OPTS)) { opts = io::configure_unsaturated_equilibrium_options( simul_data.mesh1d->nb_nodes(), simul_data.units, conf_chem.get_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_SS_EQUILIBRIUM_OPTS) ); } else { SPC_CONF_LOG << " Use default options "; opts = EquilibriumOptionsVector::make(simul_data.mesh1d->nb_nodes()); SPC_CONF_LOG << "Default SpecMiCP options : "; std::ostringstream msg; io::print_specmicp_options(msg, opts->get("default")); SPC_CONF_LOG << msg.str() << "\n" << SPC_CONF_LOG_HLINE; } SPC_CONF_LOG << SPC_CONF_LOG_HLINE; return EquilibriumStagger::make( simul_data.variables, simul_data.bcs, opts ); } std::shared_ptr initialize_transport_stagger( SimulationData& simul_data, io::YAMLConfigHandle& conf_staggers ) { SPC_CONF_LOG << "\n\n --> Transport stagger <-- \n" << SPC_CONF_LOG_SECTION; const UnsaturatedVariables* const variables = simul_data.variables.get(); const database::DataContainer* const raw_data = simul_data.raw_data.get(); // No conf // ------- if (not conf_staggers.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER)) { return UnsaturatedTransportStagger::make( simul_data.variables, simul_data.bcs, true, true); } // Conf is given // ------------- auto transport_conf = conf_staggers.get_section( SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER); const bool merge_sat = transport_conf.get_optional_attribute( SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_SATURATION, true); if (merge_sat) { SPC_CONF_LOG << "Use merged saturation/vapor pressure equation"; } else { SPC_CONF_LOG << "Use separate saturation/vapor pressure equations"; } const bool merge_aq = transport_conf.get_optional_attribute( SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_AQUEOUS, true); if (merge_aq) { SPC_CONF_LOG << "Use merged aqueous liquid/pressure equation"; } else { SPC_CONF_LOG << "Use separate aqueous liquid/pressure equations"; } auto stagger = UnsaturatedTransportStagger::make( simul_data.variables, simul_data.bcs, merge_sat, merge_aq); // ! scaling => done in variables // default options // --------------- dfpmsolver::ParabolicDriverOptions default_opts; bool rewrite_default = false; if (transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS)) { io::configure_transport_options( default_opts, transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS) ); rewrite_default = true; } SPC_CONF_LOG << "\nDefault transport options\n" << SPC_CONF_LOG_HLINE; std::ostringstream msg; io::print_transport_options(default_opts, msg); SPC_CONF_LOG << msg.str() << SPC_CONF_LOG_HLINE; // Saturation option // ----------------- auto& sat_opts = *(stagger->get_saturation_options()); sat_opts = default_opts; if (transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_SATURATION_OPTS)) { io::configure_transport_options( sat_opts, transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_SATURATION_OPTS) ); } SPC_CONF_LOG << "\nTransport options for saturation equations \n" << SPC_CONF_LOG_HLINE; msg = std::ostringstream(); io::print_transport_options(sat_opts, msg); SPC_CONF_LOG << msg.str() << SPC_CONF_LOG_HLINE; // Aqueous options // --------------- if (transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_AQUEOUS_OPTS)) { auto conf = transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_AQUEOUS_OPTS); if (not conf.is_sequence()) { // config for every aqueous equations dfpmsolver::ParabolicDriverOptions opts; if (rewrite_default) opts = default_opts; io::configure_transport_options( opts, std::move(conf)); for (index_t aq: raw_data->range_aqueous_component()) { *(stagger->get_aqueous_options(aq)) = opts; } SPC_CONF_LOG << "\nTransport options for aqueous equations \n" << SPC_CONF_LOG_HLINE; msg = std::ostringstream(); io::print_transport_options(opts, msg); SPC_CONF_LOG << msg.str() << SPC_CONF_LOG_HLINE; } else { // First set default is needed if (rewrite_default) { for (index_t aq: simul_data.raw_data->range_aqueous_component()) { *(stagger->get_aqueous_options(aq)) = default_opts; } } uindex_t size = conf.size(); for (uindex_t ind=0; ind( SPC_CF_A_COMPONENT); auto comp_id = raw_data->get_id_component(comp_label); io::configure_transport_options( *(stagger->get_aqueous_options(comp_id)), std::move(subconf) ); } } } // Gas options // ----------- if ( transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_GAS_OPTS) and not (merge_aq and merge_sat)) // skip conf if not needed { auto conf = transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_GAS_OPTS); if (not conf.is_sequence()) { // config for every gas equations dfpmsolver::ParabolicDriverOptions opts; if (rewrite_default) opts = default_opts; io::configure_transport_options( opts, std::move(conf)); if ( not merge_sat and variables->component_has_gas(0)) { *(stagger->get_gas_options(0)) = opts; } if ( not merge_aq) { for (index_t aq: raw_data->range_aqueous_component()) { if (simul_data.variables->component_has_gas(aq) ) *(stagger->get_aqueous_options(aq)) = opts; } } SPC_CONF_LOG << "\nTransport options for gaseous equations \n" << SPC_CONF_LOG_HLINE; msg = std::ostringstream(); io::print_transport_options(opts, msg); SPC_CONF_LOG << msg.str() << SPC_CONF_LOG_HLINE; } else { // First rewrite default if (rewrite_default) { // H2O if (not merge_sat and variables->component_has_gas(0)) { *(stagger->get_gas_options(0)) = default_opts; } // aqueous components if (not merge_aq) { for (index_t aq: raw_data->range_aqueous_component()) { if (simul_data.variables->component_has_gas(aq)) *(stagger->get_aqueous_options(aq)) = default_opts; } } } // Then go over all conf uindex_t size = conf.size(); for (uindex_t ind=0; ind( SPC_CF_A_COMPONENT); auto comp_id = raw_data->get_id_component(comp_label); // H2O if (comp_id == 0) { if (merge_sat) continue; io::configure_transport_options( *(stagger->get_gas_options(0)), std::move(subconf) ); } // aqueous component else { if (merge_aq) continue; io::configure_transport_options( *(stagger->get_gas_options(comp_id)), std::move(subconf) ); } } } } return stagger; } std::string initialize_output( SimulationData& simul_data, const std::string& work_dir, io::YAMLConfigHandle&& conf ) { auto type = conf.get_optional_attribute( SPC_CF_S_REACTOUTPUT_A_TYPE, SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5); if (type == SPC_CF_V_DEFAULT) type = SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5; if (type != SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5) { conf.report_error( io::YAMLConfigError::InvalidArgument, "Only hdf5 is accepted for attribute " SPC_CF_S_REACTOUTPUT_A_TYPE " for now."); } auto filepath = conf.get_required_attribute( SPC_CF_S_REACTOUTPUT_A_FILEPATH); if (not work_dir.empty() and not utils::is_path_absolute(filepath)) { filepath = utils::complete_path(work_dir, filepath); } if (conf.has_attribute(SPC_CF_S_REACTOUTPUT_A_DATABASE)) { database::Database db_manager(simul_data.raw_data); auto db_path = conf.get_attribute(SPC_CF_S_REACTOUTPUT_A_DATABASE); if (not work_dir.empty() and not utils::is_path_absolute(db_path)) { db_path = utils::complete_path(work_dir, db_path); } db_manager.save(db_path); } return filepath; } void print_db_to_conf_log(database::DataContainer* raw_data) { SPC_CONF_LOG << "Freeze database\n" << SPC_CONF_LOG_SECTION; SPC_CONF_LOG << " - " << raw_data->nb_component() << " components" << "\n - " << raw_data->nb_aqueous() << " secondary aqueous species" << "\n - " << raw_data->nb_mineral() << " solid phases at equilibrium" << "\n - " << raw_data->nb_gas() << " gas"; SPC_CONF_LOG << "List component : "; for (auto comp: raw_data->range_component()) { SPC_CONF_LOG << " - " << raw_data->get_label_component(comp); } SPC_CONF_LOG << SPC_CONF_LOG_SECTION; } diff --git a/src/reactmicp/io/configuration.cpp b/src/reactmicp/io/configuration.cpp index 1acc168..bb082d3 100644 --- a/src/reactmicp/io/configuration.cpp +++ b/src/reactmicp/io/configuration.cpp @@ -1,258 +1,260 @@ /* ============================================================================= 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 "configuration.hpp" #include "specmicp_common/io/yaml.hpp" #include "specmicp_common/io/safe_config.hpp" #include "dfpm/solver/parabolic_structs.hpp" #include "reactmicp/solver/runner.hpp" #include "reactmicp/solver/reactive_transport_solver_structs.hpp" +#include "reactmicp/solver/timestepper_structs.hpp" + #include "saturated_react.hpp" #include "specmicp_database/database.hpp" #include "specmicp_common/io/config_yaml_sections.h" #include #include namespace specmicp { namespace io { std::string clean_label(std::string label); void configure_reactmicp_options( reactmicp::solver::ReactiveTransportOptions& options, io::YAMLConfigHandle&& conf ) { conf.set_if_attribute_exists(options.residuals_tolerance, SPC_CF_S_REACTMICP_A_RES_TOL, 0.0, 1.0); conf.set_if_attribute_exists(options.absolute_residuals_tolerance, SPC_CF_S_REACTMICP_A_ABS_TOL, 0.0); conf.set_if_attribute_exists(options.step_tolerance, SPC_CF_S_REACTMICP_A_STEP_TOL, 0.0, 1.0); conf.set_if_attribute_exists(options.good_enough_tolerance, SPC_CF_S_REACTMICP_A_GOOD_ENOUGH_TOL, 0.0); conf.set_if_attribute_exists(options.implicit_upscaling, SPC_CF_S_REACTMICP_A_IMPL_UPSCALING); conf.set_if_attribute_exists(options.maximum_iterations, SPC_CF_S_REACTMICP_A_MAX_ITER, 0); } void SPECMICP_DLL_PUBLIC print_reactmicp_options( std::ostream& output, const reactmicp::solver::ReactiveTransportOptions& options ) { output << " - residuals tolerance " << options.residuals_tolerance << "\n - absolute tolerance " << options.absolute_residuals_tolerance << "\n - step tolerance " << options.step_tolerance << "\n - good enough tolerance " << options.good_enough_tolerance << "\n - maximum iterations " << options.maximum_iterations << "\n - implicit upscaling " << options.implicit_upscaling << std::endl; } void configure_reactmicp_csv_output( io::OutputNodalVariables& output_policy, const reactmicp::solver::SimulationInformation& simul_info, YAMLConfigHandle& conf ) { database::Database db_manager(output_policy.get_database()); if (conf.get_optional_attribute(SPC_CF_S_REACTOUTPUT_A_PH, false)) { output_policy.register_pH(simul_info.complete_filepath("ph", "dat")); } if (conf.get_optional_attribute(SPC_CF_S_REACTOUTPUT_A_POROSITY, false)) { output_policy.register_porosity( simul_info.complete_filepath("porosity", "dat")); } if (conf.get_optional_attribute(SPC_CF_S_REACTOUTPUT_A_DIFFUSIVITY, false)) { output_policy.register_diffusion_coefficient( simul_info.complete_filepath("diffusivity", "dat")); } if (conf.has_node(SPC_CF_S_REACTOUTPUT_A_VOLFRAC_MINERAL)) { auto labels = conf.list_to_vector( SPC_CF_S_REACTOUTPUT_A_VOLFRAC_MINERAL); for (auto label: labels) { index_t id = db_manager.safe_mineral_label_to_id(label); auto simple_label = clean_label(label); output_policy.register_volume_fraction_mineral( id, simul_info.complete_filepath("phi_"+simple_label, "dat")); } } if (conf.has_node(SPC_CF_S_REACTOUTPUT_A_TOT_CONC)) { auto labels = conf.list_to_vector( SPC_CF_S_REACTOUTPUT_A_TOT_CONC); for (auto label: labels) { index_t id = db_manager.safe_component_label_to_id(label); auto simple_label = clean_label(label); output_policy.register_total_concentration( id, simul_info.complete_filepath("t_"+simple_label, "dat")); } } if (conf.has_node(SPC_CF_S_REACTOUTPUT_A_TOT_AQ_CONC)) { auto labels = conf.list_to_vector( SPC_CF_S_REACTOUTPUT_A_TOT_AQ_CONC); for (auto label: labels) { index_t id = db_manager.safe_component_label_to_id(label); auto simple_label = clean_label(label); output_policy.register_total_aqueous_concentration( id, simul_info.complete_filepath("c_"+simple_label, "dat")); } } if (conf.has_node(SPC_CF_S_REACTOUTPUT_A_TOT_S_CONC)) { auto labels = conf.list_to_vector( SPC_CF_S_REACTOUTPUT_A_TOT_S_CONC); for (auto label: labels) { index_t id = db_manager.safe_component_label_to_id(label); auto simple_label = clean_label(label); output_policy.register_total_solid_concentration( id, simul_info.complete_filepath("s_"+simple_label, "dat")); } } if (conf.has_node(SPC_CF_S_REACTOUTPUT_A_SATINDEX)) { auto labels = conf.list_to_vector( SPC_CF_S_REACTOUTPUT_A_SATINDEX); for (auto label: labels) { index_t id = db_manager.safe_mineral_kinetic_label_to_id(label); auto simple_label = clean_label(label); output_policy.register_saturation_index_mineral_kinetic( id, simul_info.complete_filepath("si_"+simple_label, "dat")); } } } reactmicp::solver::SimulationInformation SPECMICP_DLL_PUBLIC configure_simulation_information( io::YAMLConfigHandle&& conf ) { reactmicp::solver::SimulationInformation simul_info( conf.get_required_attribute(SPC_CF_S_SIMULINFO_A_NAME), conf.get_required_attribute(SPC_CF_S_SIMULINFO_A_OUTPUTSTEP) ); conf.set_if_attribute_exists(simul_info.output_prefix, SPC_CF_S_SIMULINFO_A_OUTPUTPREFIX); conf.set_if_attribute_exists(simul_info.print_iter_info, SPC_CF_S_SIMULINFO_A_PRINTITER); return simul_info; } std::string clean_label(std::string label) { static std::vector orig {']', '[', '(', ')', '+', '-', ','}; static std::vector repl {'_', '_', '_', '_', 'p', 'm', '_'}; for (auto it = label.begin(); it!=label.end(); ++it) { auto itf = std::find(orig.begin(), orig.end(), *it); if (itf != orig.end()) { *it = repl[itf - orig.begin()]; } } std::transform(label.begin(), label.end(), label.begin(), ::tolower); return label; } void configure_reactmicp_timestepper( reactmicp::solver::TimestepperOptions& options, io::YAMLConfigHandle&& conf ) { options.lower_bound = conf.get_required_attribute( SPC_CF_S_TIMESTEP_A_LOWER_BOUND, 0.0); options.upper_bound = conf.get_required_attribute( SPC_CF_S_TIMESTEP_A_UPPER_BOUND, 0.0); if (conf.has_attribute(SPC_CF_S_TIMESTEP_A_RESTART_DT)) { options.restart_timestep = conf.get_attribute( SPC_CF_S_TIMESTEP_A_RESTART_DT, 0.0); } else { options.restart_timestep = options.lower_bound; } conf.set_if_attribute_exists(options.iteration_lower_target, SPC_CF_S_TIMESTEP_A_LOWER_TARGET, 0.0); conf.set_if_attribute_exists(options.iteration_upper_target, SPC_CF_S_TIMESTEP_A_UPPER_TARGET, 0.0); if (options.iteration_upper_target < options.iteration_lower_target) { conf.report_error( YAMLConfigError::InvalidArgument, "Iterations lower target is greater than the upper target."); } conf.set_if_attribute_exists(options.alpha_average, SPC_CF_S_TIMESTEP_A_AVERAGE_PARAMETER, 0.0, 1.0); conf.set_if_attribute_exists(options.decrease_factor, SPC_CF_S_TIMESTEP_A_FACTOR_IF_DECREASE, 0.0, 1.0); conf.set_if_attribute_exists(options.increase_factor, SPC_CF_S_TIMESTEP_A_FACTOR_IF_INCREASE, 1.0); conf.set_if_attribute_exists(options.increase_error_minimization, SPC_CF_S_TIMESTEP_A_FACTOR_IF_MINIMUM, 0.0); conf.set_if_attribute_exists(options.decrease_failure, SPC_CF_S_TIMESTEP_A_FACTOR_IF_FAILURE, 0.0); } } //end namespace io } //end namespace specmicp diff --git a/src/reactmicp/reactmicp_core.hpp b/src/reactmicp/reactmicp_core.hpp index 97468f9..6ed288d 100644 --- a/src/reactmicp/reactmicp_core.hpp +++ b/src/reactmicp/reactmicp_core.hpp @@ -1,75 +1,76 @@ /* ============================================================================= 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. * ============================================================================= */ //! \file reactmicp_core.hpp //! \brief Core files of the reactmicp solver - not standalone //! //! This file includes the main headers from the ReactMiCP solver. #ifndef SPECMICP_REACTMICP_CORE_HPP #define SPECMICP_REACTMICP_CORE_HPP #include "specmicp_common/types.hpp" #include "specmicp_database/database.hpp" #include "specmicp/specmicp.hpp" #include "dfpm/mesh.hpp" // solver #include "reactmicp/solver/reactive_transport_solver.hpp" #include "reactmicp/solver/reactive_transport_solver_structs.hpp" #include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" #include "reactmicp/solver/timestepper.hpp" +#include "reactmicp/solver/timestepper_structs.hpp" #include "reactmicp/solver/runner.hpp" // output #include "specmicp_common/io/csv_formatter.hpp" #include "dfpm/io/meshes.hpp" #include "reactmicp/io/reactive_transport.hpp" #include "specmicp_common/physics/io/units.hpp" // configuration #include "specmicp_common/io/yaml.hpp" #include "specmicp_database/io/configuration.hpp" #include "dfpm/io/configuration.hpp" #include "reactmicp/io/configuration.hpp" #include "specmicp_common/timer.hpp" #endif // SPECMICP_REACTMICP_CORE_HPP diff --git a/src/reactmicp/solver/CMakeLists.txt b/src/reactmicp/solver/CMakeLists.txt index c6bdf0e..fd115cf 100644 --- a/src/reactmicp/solver/CMakeLists.txt +++ b/src/reactmicp/solver/CMakeLists.txt @@ -1,22 +1,23 @@ set( reactmicp_solver_srcs reactive_transport_solver.cpp runner.cpp timestepper.cpp ) set( reactmicp_solver_headers reactive_transport_solver.hpp reactive_transport_solver_structs.hpp runner.hpp timestepper.hpp + timestepper_structs.hpp ) add_to_reactmicp_srcs_list(reactmicp_solver_srcs) add_to_reactmicp_headers_list(reactmicp_solver_headers) INSTALL(FILES ${reactmicp_solver_headers} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/solver ) add_subdirectory( staggers_base ) diff --git a/src/reactmicp/solver/runner.cpp b/src/reactmicp/solver/runner.cpp index 35aafb3..67e52ad 100644 --- a/src/reactmicp/solver/runner.cpp +++ b/src/reactmicp/solver/runner.cpp @@ -1,149 +1,228 @@ /* ============================================================================= 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 "runner.hpp" +#include "timestepper.hpp" +#include "timestepper_structs.hpp" +#include "reactive_transport_solver.hpp" +#include "reactive_transport_solver_structs.hpp" + #include "reactmicp/io/reactive_transport.hpp" #include "staggers_base/variables_base.hpp" #include "specmicp_common/timer.hpp" #include "specmicp_common/io/csv_formatter.hpp" #include "specmicp_common/filesystem.hpp" #include namespace specmicp { namespace reactmicp { namespace solver { -void ReactiveTransportRunner::run_until(scalar_t target, VariablesBasePtr variables) +struct ReactiveTransportRunner::ReactiveTransportRunnerImpl +{ + ReactiveTransportRunnerImpl( + ReactiveTransportSolver& solver, + scalar_t lower_dt_bound, + scalar_t upper_dt_bound, + const SimulationInformation& info + ): + m_solver(solver), + m_timestepper(lower_dt_bound, upper_dt_bound, 0, 2.0), + m_simulinfo(info) + { + + } + + index_t m_cnt{0}; + ReactiveTransportSolver& m_solver; + Timestepper m_timestepper; + const SimulationInformation& m_simulinfo; + output_f m_output_function{dummy_output}; + scalar_t m_output_target; + + scalar_t run_until(scalar_t target, VariablesBasePtr variables); +}; + +ReactiveTransportRunner::ReactiveTransportRunner( + ReactiveTransportSolver& solver, + scalar_t lower_dt_bound, + scalar_t upper_dt_bound, + const SimulationInformation& info): + m_impl(utils::make_pimpl( + solver, lower_dt_bound, upper_dt_bound, info) + ) +{} + + +ReactiveTransportRunner::~ReactiveTransportRunner() = default; + + +void ReactiveTransportRunner::set_output_policy(output_f output_policy) { + m_impl->m_output_function = output_policy; +} + +ReactiveTransportOptions& ReactiveTransportRunner::get_options() { + return m_impl->m_solver.get_options(); +} + +TimestepperOptions& ReactiveTransportRunner::get_timestepper_options() { + return m_impl->m_timestepper.get_options(); +} + +ReactiveTransportPerformance& ReactiveTransportRunner::get_perfs() { + return m_impl->m_solver.get_perfs(); +} + +ReactiveTransportTimer& ReactiveTransportRunner::get_timer() { + return m_impl->m_solver.get_timer(); +} + +scalar_t ReactiveTransportRunner::run_until( + scalar_t target, + VariablesBasePtr variables + ) +{ + return m_impl->run_until(target, variables); +} + + +// Implementation + +scalar_t ReactiveTransportRunner::ReactiveTransportRunnerImpl::run_until( + scalar_t target, + VariablesBasePtr variables + ) { Timer total_time; VariablesBase* vars = variables.get(); total_time.start(); m_timestepper.set_total_target(target); m_output_target = m_timestepper.get_total()+m_simulinfo.output_step; auto filepath = m_simulinfo.complete_filepath( std::string("iter"), std::string("dat")); if (not m_simulinfo.working_dir.empty() and not utils::is_path_absolute(filepath)) { filepath = utils::complete_path(m_simulinfo.working_dir, filepath); } io::OutFile out_iter(filepath); io::print_reactmicp_header(out_iter); if (m_simulinfo.print_iter_info) { io::print_reactmicp_performance_long_header(out_iter); } scalar_t dt = m_timestepper.get_options().restart_timestep; auto retcode = ReactiveTransportReturnCode::NotConvergedYet; bool last_solved = true; while(retcode >= ReactiveTransportReturnCode::NotConvergedYet and m_timestepper.get_total() < m_timestepper.get_total_target()) { if (m_timestepper.get_total() > 0) last_solved = false; Timer step_timer; step_timer.start(); retcode = m_solver.solve_timestep(dt, vars); step_timer.stop(); if (m_simulinfo.print_iter_info) io::print_reactmicp_performance_long( out_iter, m_cnt, m_timestepper.get_total()+dt, m_solver.get_perfs() ); dt = m_timestepper.next_timestep( dt, retcode, m_solver.get_perfs().nb_iterations ); if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) { dt = m_timestepper.get_options().restart_timestep; vars->reset_main_variables(); retcode = m_solver.solve_timestep(dt, vars); if (m_simulinfo.print_iter_info) io::print_reactmicp_performance_long( out_iter, m_cnt, m_timestepper.get_total()+dt, m_solver.get_perfs() ); if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) vars->reset_main_variables(); dt = m_timestepper.next_timestep( dt, retcode, m_solver.get_perfs().nb_iterations ); } ++m_cnt; // output if (m_timestepper.get_total() > m_output_target) { last_solved = true; m_output_function(m_timestepper.get_total(), variables); m_output_target += m_simulinfo.output_step; } if (retcode == reactmicp::solver::ReactiveTransportReturnCode::UserTermination) break; } if (not last_solved) { m_output_function(m_timestepper.get_total(), variables); } total_time.stop(); - io::print_reactmicp_end(out_iter, total_time, get_timer()); + io::print_reactmicp_end(out_iter, total_time, m_solver.get_timer()); + return m_timestepper.get_total(); } } // end namespace solver } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/solver/runner.hpp b/src/reactmicp/solver/runner.hpp index 0d6435e..8a27bb9 100644 --- a/src/reactmicp/solver/runner.hpp +++ b/src/reactmicp/solver/runner.hpp @@ -1,134 +1,146 @@ /* ============================================================================= 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_REACTMICP_SOLVER_RUNNER_HPP #define SPECMICP_REACTMICP_SOLVER_RUNNER_HPP -#include "reactive_transport_solver.hpp" -#include "timestepper.hpp" +//! \file reactmicp/solver/runner.hpp +//! \brief Run the reactive transport solver until a target is reached + +#include "specmicp_common/types.hpp" +#include "specmicp_common/pimpl_ptr.hpp" + #include #include +#include namespace specmicp { namespace reactmicp { namespace solver { +// forward declarations +class ReactiveTransportSolver; +struct ReactiveTransportOptions; +struct ReactiveTransportPerformance; +struct ReactiveTransportTimer; +struct TimestepperOptions; + +class VariablesBase; +using VariablesBasePtr = std::shared_ptr; //! \struct SimulationInformation //! \brief Information about the simulation struct SPECMICP_DLL_PUBLIC SimulationInformation { std::string name; //!< Name of the simulation std::string output_prefix; //!< prefix for the output files std::string working_dir {""}; bool print_iter_info{true}; //!< If true, print the iteration informations scalar_t output_step; //!< output step SimulationInformation(std::string name_simul, scalar_t outputstep): name(name_simul), output_prefix(name_simul+"_"), output_step(outputstep) - {} - std::string complete_filepath(const std::string& name, const std::string& suffix) const { + //! \brief Complete the name of an output file + std::string complete_filepath( + const std::string& name, + const std::string& suffix + ) const { return output_prefix+name+"."+suffix; } - std::string complete_filepath(std::string&& name, std::string&& suffix) const { + //! \brief Complete the name of an output file + std::string complete_filepath( + std::string&& name, + std::string&& suffix + ) const { return output_prefix+name+"."+suffix; } }; - +//! \brief Signature of an output function using output_f = std::function; - +//! \brief No output = default inline void dummy_output(scalar_t _, VariablesBasePtr __) {} +//! \brief Run the reactive transport solver until a target class SPECMICP_DLL_PUBLIC ReactiveTransportRunner { public: + //! + //! \param solver the reactive transport solver + //! \param lower_dt_bound the timestep lower bound + //! \param upper_dt_bound the timestep upper bound + //! \param info information about the simulation ReactiveTransportRunner(ReactiveTransportSolver& solver, scalar_t lower_dt_bound, scalar_t upper_dt_bound, - const SimulationInformation& info): - m_solver(solver), - m_timestepper(lower_dt_bound, upper_dt_bound, 0, 2.0), - m_simulinfo(info) - { - - } - + const SimulationInformation& info); + ~ReactiveTransportRunner(); - void run_until(scalar_t target, VariablesBasePtr variables); + //! \brief Run the solver until a certain target + //! + //! \param target the target to reach + //! \param variables the variables to use + scalar_t run_until(scalar_t target, VariablesBasePtr variables); - void set_output_policy(output_f output_policy) { - m_output_function = output_policy; - } + //! \brief Set the output function + void set_output_policy(output_f output_policy); - ReactiveTransportOptions& get_options() { - return m_solver.get_options(); - } + //! \brief The options to the reactive transport solver + ReactiveTransportOptions& get_options(); - TimestepperOptions& get_timestepper_options() { - return m_timestepper.get_options(); - } + //! \brief The options to the timestepper + TimestepperOptions& get_timestepper_options(); - ReactiveTransportPerformance& get_perfs() { - return m_solver.get_perfs(); - } + //! \brief The performance strcts of the reactive transport solver + ReactiveTransportPerformance& get_perfs(); - ReactiveTransportTimer& get_timer() { - return m_solver.get_timer(); - } + //! \brief Return the timer of the reactive transport solver + ReactiveTransportTimer& get_timer(); private: - index_t m_cnt{0}; - ReactiveTransportSolver& m_solver; - Timestepper m_timestepper; - const SimulationInformation& m_simulinfo; - - output_f m_output_function{dummy_output}; - - scalar_t m_output_target; - - + struct SPECMICP_DLL_LOCAL ReactiveTransportRunnerImpl; + utils::pimpl_ptr m_impl; }; } // end namespace solver } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SOLVER_RUNNER_HPP diff --git a/src/reactmicp/solver/timestepper.hpp b/src/reactmicp/solver/timestepper.hpp index 0cd1737..584e841 100644 --- a/src/reactmicp/solver/timestepper.hpp +++ b/src/reactmicp/solver/timestepper.hpp @@ -1,121 +1,115 @@ /* ============================================================================= 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_REACTMICP_SOLVER_TIMESTEPPER_HPP #define SPECMICP_REACTMICP_SOLVER_TIMESTEPPER_HPP +//! \file reactmicp/solver/timestepper.hpp +//! \brief Manage the timestep duration in the reactive transport solver + #include "specmicp_common/moving_average.hpp" #include "specmicp_common/options_handler.hpp" +#include "timestepper_structs.hpp" + namespace specmicp { namespace reactmicp { namespace solver { // forward declaration enum class ReactiveTransportReturnCode; -//! \brief Options for the timestepper -struct SPECMICP_DLL_PUBLIC TimestepperOptions -{ - scalar_t lower_bound; //!< Lower bound for the timestep - scalar_t upper_bound; //!< Upper bound for the timestep - scalar_t restart_timestep; //!< Value used when restarting the problem - - scalar_t iteration_lower_target{1.01}; //!< Lower target for the number of iterations - scalar_t iteration_upper_target{15.0}; //!< Upper target for the number of iterations - - scalar_t alpha_average{0.5}; //!< Parameter for the exponential moving average - - scalar_t decrease_failure{0.5}; //!< Reduction factor in case of failure - scalar_t increase_error_minimization{1.3}; //!< Increase factor in case of error minimization - scalar_t decrease_factor{0.75}; //!< Reduction factor to get the number of iterations inside the target - scalar_t increase_factor{1.25}; //!< Increase factor to get the number of iterations inside the target - - TimestepperOptions(scalar_t dt_lower_bound, scalar_t dt_upper_bound): - lower_bound(dt_lower_bound), - upper_bound(dt_upper_bound), - restart_timestep(dt_lower_bound) - {} -}; //! \brief Adaptative timestepper for the reactive transport solver class SPECMICP_DLL_PUBLIC Timestepper: public OptionsHandler { public: //! \brief Constructor //! //! \param dt_lower_bound lower_bound for the timestep //! \param dt_upper_bound upper_bound for the timestep //! \param total_target total target time //! \param init_iterations initial iterations Timestepper(scalar_t dt_lower_bound, scalar_t dt_upper_bound, scalar_t total_target, scalar_t init_iterations ): OptionsHandler(dt_lower_bound, dt_upper_bound), m_total(0), m_total_target(total_target), m_average(0.5, init_iterations) { m_average.set_alpha(get_options().alpha_average); } //! \brief Return the total time scalar_t get_total() const {return m_total;} //! \brief Return the total target time scalar_t get_total_target() const {return m_total_target;} //! \brief Set the total target time void set_total_target(scalar_t total_target) {m_total_target = total_target;} + //! \brief Set the average parameter + //! + //! \param alpha parameter of the exponential moving average void set_average_parameter(scalar_t alpha) { get_options().alpha_average = alpha; m_average.set_alpha(alpha); } - //! obtain the next timestep - scalar_t next_timestep(scalar_t dt_done, ReactiveTransportReturnCode return_code, index_t nb_iterations); + //! \brief obtain the next timestep + //! + //! \param dt_done duration of the previous timestep + //! \param return_code return code of the previous timestep + //! \param nb_iterations number of fix-point iterations of the previous timestep + //! \return a timestep duration + scalar_t next_timestep( + scalar_t dt_done, + ReactiveTransportReturnCode return_code, + index_t nb_iterations + ); private: scalar_t m_total; scalar_t m_total_target; utils::ExponentialMovingAverage m_average; }; } // end namespace solver } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SOLVER_TIMESTEPPER_HPP diff --git a/src/reactmicp/solver/timestepper.hpp b/src/reactmicp/solver/timestepper_structs.hpp similarity index 58% copy from src/reactmicp/solver/timestepper.hpp copy to src/reactmicp/solver/timestepper_structs.hpp index 0cd1737..a520046 100644 --- a/src/reactmicp/solver/timestepper.hpp +++ b/src/reactmicp/solver/timestepper_structs.hpp @@ -1,121 +1,75 @@ /* ============================================================================= 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_REACTMICP_SOLVER_TIMESTEPPER_HPP -#define SPECMICP_REACTMICP_SOLVER_TIMESTEPPER_HPP +#ifndef SPECMICP_REACTMICP_SOLVER_TIMESTEPPER_STRUCTS_HPP +#define SPECMICP_REACTMICP_SOLVER_TIMESTEPPER_STRUCTS_HPP -#include "specmicp_common/moving_average.hpp" -#include "specmicp_common/options_handler.hpp" +//! \file reactmicp/solver/timestepper_structs +//! \brief Options for the timestepper + +#include "specmicp_common/types.hpp" namespace specmicp { namespace reactmicp { namespace solver { -// forward declaration -enum class ReactiveTransportReturnCode; - //! \brief Options for the timestepper struct SPECMICP_DLL_PUBLIC TimestepperOptions { - scalar_t lower_bound; //!< Lower bound for the timestep - scalar_t upper_bound; //!< Upper bound for the timestep + scalar_t lower_bound; //!< Lower bound for the timestep + scalar_t upper_bound; //!< Upper bound for the timestep scalar_t restart_timestep; //!< Value used when restarting the problem scalar_t iteration_lower_target{1.01}; //!< Lower target for the number of iterations scalar_t iteration_upper_target{15.0}; //!< Upper target for the number of iterations scalar_t alpha_average{0.5}; //!< Parameter for the exponential moving average - scalar_t decrease_failure{0.5}; //!< Reduction factor in case of failure + scalar_t decrease_failure{0.5}; //!< Reduction factor in case of failure scalar_t increase_error_minimization{1.3}; //!< Increase factor in case of error minimization scalar_t decrease_factor{0.75}; //!< Reduction factor to get the number of iterations inside the target scalar_t increase_factor{1.25}; //!< Increase factor to get the number of iterations inside the target TimestepperOptions(scalar_t dt_lower_bound, scalar_t dt_upper_bound): lower_bound(dt_lower_bound), upper_bound(dt_upper_bound), restart_timestep(dt_lower_bound) {} }; -//! \brief Adaptative timestepper for the reactive transport solver -class SPECMICP_DLL_PUBLIC Timestepper: public OptionsHandler -{ -public: - //! \brief Constructor - //! - //! \param dt_lower_bound lower_bound for the timestep - //! \param dt_upper_bound upper_bound for the timestep - //! \param total_target total target time - //! \param init_iterations initial iterations - Timestepper(scalar_t dt_lower_bound, - scalar_t dt_upper_bound, - scalar_t total_target, - scalar_t init_iterations - ): - OptionsHandler(dt_lower_bound, dt_upper_bound), - m_total(0), - m_total_target(total_target), - m_average(0.5, init_iterations) - { - m_average.set_alpha(get_options().alpha_average); - } - - //! \brief Return the total time - scalar_t get_total() const {return m_total;} - //! \brief Return the total target time - scalar_t get_total_target() const {return m_total_target;} - //! \brief Set the total target time - void set_total_target(scalar_t total_target) {m_total_target = total_target;} - - void set_average_parameter(scalar_t alpha) { - get_options().alpha_average = alpha; - m_average.set_alpha(alpha); - } - - //! obtain the next timestep - scalar_t next_timestep(scalar_t dt_done, ReactiveTransportReturnCode return_code, index_t nb_iterations); - -private: - scalar_t m_total; - scalar_t m_total_target; - - utils::ExponentialMovingAverage m_average; -}; - } // end namespace solver } // end namespace reactmicp } // end namespace specmicp -#endif // SPECMICP_REACTMICP_SOLVER_TIMESTEPPER_HPP +#endif // SPECMICP_REACTMICP_SOLVER_TIMESTEPPER_STRUCTS_HPP