diff --git a/examples/reactmicp/saturated_react/carbonation.cpp b/examples/reactmicp/saturated_react/carbonation.cpp index d7dde40..94f6748 100644 --- a/examples/reactmicp/saturated_react/carbonation.cpp +++ b/examples/reactmicp/saturated_react/carbonation.cpp @@ -1,576 +1,576 @@ /* ============================================================================= 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. * ============================================================================= */ // ====================================================== // Leaching of a cement paste in CO2-saturated water // ======================================================= // This file can be used to learn how to use ReactMiCP // Include ReactMiCP #include "reactmicp/reactmicp.hpp" // Other includes that would be needed later : #include #include "specmicp_common/physics/laws.hpp" #include "specmicp_common/timer.hpp" #include "specmicp_common/io/safe_config.hpp" #include "specmicp_common/physics/io/configuration.hpp" #include "specmicp_common/io/config_yaml_sections.h" #include #include #include // we use the following namespaces using namespace specmicp; using namespace reactmicp; using namespace reactmicp::systems; using namespace specmicp::reactmicp::systems::satdiff; using namespace reactmicp::solver; // We declare some functions that we will be using : // Boundary and initial conditions // =============================== // Return the chemistry constraints for the inital condition AdimensionalSystemConstraints get_cement_initial_constraints( RawDatabasePtr the_database, const units::UnitsSet& units_set); // return the chemistry constraints for the boundary condition AdimensionalSystemConstraints get_bc_initial_constraints( RawDatabasePtr the_database, const units::UnitsSet& units_set); Vector initial_compo(const Vector& publi); void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species); // The mesh // ======== // There is two version // a uniform mesh (cf config) // and a non uniform mesh mesh::Mesh1DPtr get_mesh(); // =================== // MAIN // =================== int main() { // initialize eigen // This is needed if OpenMP is used Eigen::initParallel(); // Setup the log // The output will be on cerr specmicp::logger::ErrFile::stream() = &std::cerr; // And only Warning and error messages will be printed specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; // Configuration // auto config = specmicp::io::parse_yaml_file("carbonation.yaml"); auto config = specmicp::io::YAMLConfigFile::make("carbonation.yaml"); // The database // ------------- // Obtain the database RawDatabasePtr the_database = specmicp::io::configure_database( config->get_section(SPC_CF_S_DATABASE) ); // SpecMiCP options // ---------------- // Set the options AdimensionalSystemSolverOptions specmicp_options; units::UnitsSet units_set = specmicp::io::configure_units( config->get_section(SPC_CF_S_UNITS)); specmicp::io::configure_specmicp_options(specmicp_options, units_set, config->get_section(SPC_CF_S_SPECMICP)); // The mesh // -------- // Obtain the mesh mesh::Mesh1DPtr the_mesh = specmicp::io::configure_mesh( config->get_section(SPC_CF_S_MESH)); index_t nb_nodes = the_mesh->nb_nodes(); // just boilerplate to make things easier // Print the mesh io::print_mesh("out_carbo_mesh.dat", the_mesh, units::LengthUnit::centimeter); // Initial and boundary conditions // -------------------------------- // Obtain the initial constraints // Note : the database is modified at this point // Only the relevant components are kept at this point AdimensionalSystemConstraints cement_constraints = get_cement_initial_constraints(the_database, units_set); // Freeze the database the_database->freeze_db(); // Obtain the initial condition AdimensionalSystemSolution cement_solution = solve_equilibrium( the_database, cement_constraints, specmicp_options); if(not the_database->is_valid()) throw std::runtime_error("The database has been changed during the initial condition computation"); database::Database db_manager(the_database); db_manager.save("out_carbo_db.js"); // Just to check that everything is ok ! AdimensionalSystemSolutionExtractor extractor(cement_solution, the_database, units_set); std::cout << "Porosity : " << extractor.porosity() << std::endl; std::cout << "Water volume fraction " << extractor.volume_fraction_water() << std::endl; //return EXIT_SUCCESS; // Obtain the boundary conditions AdimensionalSystemConstraints bc_constraints = get_bc_initial_constraints( the_database, specmicp_options.units_set ); AdimensionalSystemSolution bc_solution = solve_equilibrium( the_database, bc_constraints, specmicp_options ); if(not the_database->is_valid()) throw std::runtime_error("The database has been changed during the boundary conditions computation"); // The constraints // -------------- // Boundary condition // 'list_fixed_nodes' lists the node with fixed composition std::vector list_fixed_nodes = {0}; // list_constraints lists the different speciation constraints in the system // Note : the total concentrations are computed from the initial solution, // and thus different total concentrations do not constitute a reason to create different constraints std::vector list_constraints = {cement_constraints, bc_constraints}; // index_constraints links a node to the set of constraints that will be used std::vector index_constraints(nb_nodes, 0); index_constraints[0] = 1; // Boundary conditions // This is the list of initial composition std::vector list_initial_composition = {cement_solution, bc_solution}; // This list links a node to its initial conditions std::vector index_initial_state(nb_nodes, 0); index_initial_state[0] = 1; // boundary condition // Variables // ----------- // Create the main set of variables // They will be initialized using the list of initial composition satdiff::SaturatedVariablesPtr variables = satdiff::init_variables(the_mesh, the_database, units_set, list_fixed_nodes, list_initial_composition, index_initial_state); // Staggers // --------- // This section initialize the staggers // The equilibrium stagger // - - - - - - - - - - - - std::shared_ptr equilibrium_stagger = std::make_shared( list_constraints, index_constraints, specmicp_options); // The transport stagger // - - - - - - - - - - - std::vector list_fixed_component {0, 1, the_database->get_id_component("Si(OH)4")}; std::map slave_nodes {}; std::shared_ptr transport_stagger = std::make_shared( variables, list_fixed_nodes, slave_nodes, list_fixed_component ); auto staggers = config->get_section(SPC_CF_S_STAGGERS); auto transport_stagger_conf = staggers.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER); specmicp::io::configure_transport_options( transport_stagger->options_solver(), transport_stagger_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS)); // The upscaling stagger // - - - - - - - - - - - CappedPowerLawParameters upscaling_param; upscaling_param.d_eff_0 = 2.726e-12*1e4; upscaling_param.cap = 1e10*1e6; upscaling_param.exponent = 3.32; upscaling_param.porosity_0 = 0.25; upscaling_param.porosity_res = 0.02; CappedPowerLaw upscaling_law(upscaling_param); UpscalingStaggerPtr upscaling_stagger = std::make_shared(upscaling_law.get_law(),the_database, units_set); // Solver // ------ // This section initialize the reactive transport solver ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger); specmicp::io::configure_reactmicp_options(solver.get_options(), config->get_section(SPC_CF_S_REACTMICP)); transport_stagger->initialize(variables); equilibrium_stagger->initialize(variables); upscaling_stagger->initialize(variables); // Some basic information about the simulation SimulationInformation simul_info = specmicp::io::configure_simulation_information( config->get_section(SPC_CF_S_SIMULINFO)); // Runner // ------- // Create the runner : loop through the timestep ReactiveTransportRunner runner(solver, 0.0, 0.0, simul_info); specmicp::io::configure_reactmicp_timestepper( runner.get_timestepper_options(), config->get_section(SPC_CF_S_TIMESTEPPER) ); // Output // ------ io::OutputNodalVariables output_policy(the_database, the_mesh, specmicp_options.units_set); if (config->has_section(SPC_CF_S_REACTOUTPUT)) { auto conf = config->get_section(SPC_CF_S_REACTOUTPUT); io::configure_reactmicp_csv_output(output_policy, simul_info, conf); } // and bind it to the runner runner.set_output_policy(output_policy.get_output_for_reactmicp()); // Solve the problem // ----------------- auto run_config = config->get_section(SPC_CF_S_RUN); - auto run_until = run_config.get_required_attribute(SPC_CF_S_RUN_A_RUNUTNIL); + auto run_until = run_config.get_required_attribute(SPC_CF_S_RUN_A_RUNUNTIL); runner.run_until(run_until, cast_var_to_base(variables)); if(not the_database->is_valid()) throw std::runtime_error("The database has been changed during the computation"); io::AdimensionalSystemSolutionSaver saver(the_database); saver.save_solutions(simul_info.complete_filepath("solutions", "dat"), variables->equilibrium_solutions()); // output information at the end of the timestep io::print_minerals_profile(the_database, variables, the_mesh, specmicp_options.units_set, simul_info.complete_filepath("profiles_end", "dat")); io::print_components_total_aqueous_concentration(the_database, variables, the_mesh, specmicp_options.units_set, simul_info.complete_filepath("c_profiles_end", "dat")); io::print_components_total_solid_concentration(the_database, variables, the_mesh, specmicp_options.units_set, simul_info.complete_filepath("s_profiles_end", "dat")); return EXIT_SUCCESS; } AdimensionalSystemConstraints get_cement_initial_constraints( RawDatabasePtr the_database, const units::UnitsSet& units_set) { Vector oxyde_compo(4); oxyde_compo << 66.98, 21.66, 7.19, 2.96; Vector species_compo; compo_from_oxyde(oxyde_compo, species_compo); scalar_t mult = 1e-6; //species_compo *= 0.947e-2; // poro 0.47 //species_compo *= 1.08e-2; // poro 0.4 //species_compo *= 1.25e-2; // poro 0.3 species_compo *= 1.1e-2; Formulation formulation; formulation.mass_solution = 1.0; formulation.amount_minerals = { {"C3S", species_compo(0)}, {"C2S", species_compo(1)}, {"C3A", species_compo(2)}, {"Gypsum", species_compo(3)}, {"Calcite", species_compo(0)*1e-3} }; formulation.extra_components_to_keep = {"HCO3[-]"}; formulation.concentration_aqueous = { {"NaOH", mult*1.0298}, {"KOH", mult*0.0801}, {"Cl[-]", mult*0.0001} }; Dissolver dissolver(the_database); dissolver.set_units(units_set); AdimensionalSystemConstraints constraints; constraints.set_charge_keeper(the_database->get_id_component("HO[-]")); constraints.total_concentrations = dissolver.dissolve(formulation); std::cout << constraints.total_concentrations << std::endl; constraints.set_saturated_system(); return constraints; } AdimensionalSystemConstraints get_bc_initial_constraints( RawDatabasePtr the_database, const units::UnitsSet& units_set) { const scalar_t rho_w = laws::density_water(units::celsius(25.0), units_set.length, units_set.mass); const scalar_t nacl = 0.3*rho_w; const scalar_t kcl = 0.28*rho_w; Vector total_concentrations; total_concentrations.setZero(the_database->nb_component()); total_concentrations(the_database->get_id_component("K[+]")) = kcl; total_concentrations(the_database->get_id_component("Cl[-]")) = nacl + kcl; total_concentrations(the_database->get_id_component("Na[+]")) = nacl; //total_concentrations(the_database->get_id_component("Si(OH)4")) = 0.0001*rho_w; AdimensionalSystemConstraints constraints; constraints.add_fixed_activity_component(the_database->get_id_component("HO[-]"), -10.3); constraints.set_charge_keeper(the_database->get_id_component("HCO3[-]")); constraints.total_concentrations = total_concentrations; constraints.set_saturated_system(); constraints.disable_surface_model(); constraints.fixed_fugacity_cs ={}; return constraints; } AdimensionalSystemSolution get_bc_initial_compo( RawDatabasePtr the_database, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ) { Vector variables; AdimensionalSystemSolver solver(the_database, constraints, options); micpsolver::MiCPPerformance perf = solver.solve(variables, true); if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet) throw std::runtime_error("Failed to solve initial composition"); std::cout << variables << std::endl; return solver.get_raw_solution(variables); } void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species) { constexpr double M_CaO = 56.08; constexpr double M_SiO2 = 60.09; constexpr double M_Al2O3 = 101.96; constexpr double M_SO3 = 80.06; Eigen::MatrixXd compo_in_oxyde(4, 4); Vector molar_mass(4); molar_mass << 1.0/M_CaO, 1.0/M_SiO2, 1.0/M_Al2O3, 1.0/M_SO3; compo_oxyde = compo_oxyde.cwiseProduct(molar_mass); //C3S C2S C3A Gypsum compo_in_oxyde << 3, 2, 3, 1, // CaO 1, 1, 0, 0, // SiO2 0, 0, 1, 0, // Al2O3 0, 0, 0, 1; // SO3 Eigen::ColPivHouseholderQR solver(compo_in_oxyde); compo_species = solver.solve(compo_oxyde); } Vector initial_compo(const Vector& publi) { Matrix publi_stochio(5,5); publi_stochio << 1, 0, 0, 0, 0, 9, 6, 0, 0, 0, 3, 0, 1, 1, 0, 3, 0, 1, 3, 0, 0, 0, 0, 0, 1; Matrix cemdata_stochio(5, 5); cemdata_stochio << 1, 0, 0, 0, 0, 1.0/0.6, 1.0, 0, 0, 0, 3, 0, 1, 1, 0, 3, 0, 1, 3, 0, 0, 0, 0, 0, 1; Vector oxyde; Eigen::ColPivHouseholderQR solver(publi_stochio); oxyde = solver.solve(publi); Vector for_cemdata = cemdata_stochio*oxyde; return for_cemdata; } mesh::Mesh1DPtr get_mesh() { Vector radius(71); const scalar_t height = 20; radius << 3.751, 3.750, 3.745, 3.740, 3.735, 3.730, 3.720, 3.710, 3.700, 3.690, 3.680, 3.670, 3.660, 3.645, 3.630, 3.615, 3.600, 3.580, 3.560, 3.540, 3.520, 3.500, 3.475, 3.450, 3.425, 3.400, 3.375, 3.350, 3.325, 3.300, 3.275, 3.250, 3.225, 3.200, 3.175, 3.150, 3.125, 3.100, 3.050, 3.025, 3.000, 2.950, 2.900, 2.850, 2.800, 2.750, 2.700, 2.650, 2.600, 2.550, 2.500, 2.450, 2.400, 2.350, 2.300, 2.250, 2.200, 2.150, 2.100, 2.050, 2.000, 1.950, 1.900, 1.850, 1.800, 1.750, 1.700, 1.650, 1.600, 1.550, 1.500 ; radius = radius/10; return mesh::axisymmetric_mesh1d(radius, height); } diff --git a/src/bin/reactmicp/unsaturated.cpp b/src/bin/reactmicp/unsaturated.cpp index 5e3ca8a..c5a0d4d 100644 --- a/src/bin/reactmicp/unsaturated.cpp +++ b/src/bin/reactmicp/unsaturated.cpp @@ -1,954 +1,954 @@ /* ============================================================================= 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 "dfpm/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" #include "specmicp_common/io/all_io_files.hpp" #define IO_FILE_OUTPUT "Main" #define IO_FILE_INIT_DB "Raw database" #define IO_FILE_WORK_DB "Work database" #define IO_ITER_FILE "Iterations log" #define IO_INPUT "IO" 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, io::AllIOFiles& all_io_files ); 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, io::AllIOFiles& all_io_file ); io::AllIOFiles initialize_all_io( 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.add_option('o', "overwrite", false, "If true, allow overwrite of previous simulation"); parser.add_option('c', "clean", false, "Clean previous output files"); parser.add_option('p', "pretend", false, "Dry-run, do not run simulation"); 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"); bool allow_overwrite = parser.get_option("overwrite"); bool clean_outputs = parser.get_option("clean"); bool only_pretend = parser.get_option("pretend"); std::vector db_dirs; if (not database_dir.empty()) {db_dirs.push_back(database_dir);} db_dirs.push_back(utils::get_current_directory()); io::add_db_dirs_from_env(db_dirs); if (not working_dir.empty()) {db_dirs.push_back(working_dir);} // Configuration File // ------------------ auto config = io::YAMLConfigFile::make(input_filepath); solver::SimulationInformation simul_info = io::configure_simulation_information( config->get_section(SPC_CF_S_SIMULINFO)); simul_info.working_dir = working_dir; auto io_file_basename = simul_info.output_prefix + "io.yml"; auto io_file_path = utils::complete_path(working_dir, io_file_basename); io::AllIOFilesMode flag = io::AllIOFilesMode::ErrorIfExist; if (allow_overwrite) { flag = io::AllIOFilesMode::Write; } else if (clean_outputs ){ flag = io::AllIOFilesMode::Read; } std::unique_ptr all_io_files = make_unique(io_file_path, flag); all_io_files->add_configuration_file(io::input_file(IO_INPUT, input_filepath)); if (clean_outputs) { if (only_pretend) { return EXIT_SUCCESS; } all_io_files->clean_output_files(); all_io_files.reset(nullptr); return EXIT_SUCCESS; // just clean, nothing to run } 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), all_io_files.get() // log file registered here ); } 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), all_io_files.get() // log file registered here ); } initialize_plugin_manager(*(config.get())); // openmp initialization // --------------------- utils::initialize_parallel(); // Simulation Data // --------------- SimulationData simul_data = initialize_simul_data( *config.get(), db_dirs, *all_io_files); all_io_files->sync(); // 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::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), *all_io_files); 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); } const auto iter_path = runner.get_iter_file_path(); if (not iter_path.empty()) { all_io_files->add_log_file( io::output_file(IO_ITER_FILE, iter_path)); } all_io_files->sync(); // 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); + SPC_CF_S_RUN_A_RUNUNTIL); // At this point we assume it's safe to destroy the config // to get information about unread keys config.reset(nullptr); if (not only_pretend) { runner.run_until(run_until, simul_data.variables); } all_io_files.reset(nullptr); return EXIT_SUCCESS; } // ==================== // // Simulation data // // ==================== // SimulationData initialize_simul_data( io::YAMLConfigHandle& safe_config, const std::vector& db_dirs, io::AllIOFiles& all_io_files ) { // 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 ); auto& metadata = simul_data.raw_data->metadata; all_io_files.add_database( io::input_file(IO_FILE_INIT_DB, metadata.path, metadata.name + " - " + metadata.version )); 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 ); } // Conf is given // ------------- UnsaturatedTransportStaggerOptions opts; 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, false); if (merge_sat) { SPC_CONF_LOG << "Use merged saturation/vapor pressure equation"; } else { SPC_CONF_LOG << "Use separate saturation/vapor pressure equations"; } opts.merge_saturation_pressure = merge_sat; const bool merge_aq = transport_conf.get_optional_attribute( SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_AQUEOUS, false); if (merge_aq) { SPC_CONF_LOG << "Use merged aqueous liquid/pressure equation"; } else { SPC_CONF_LOG << "Use separate aqueous liquid/pressure equations"; } opts.merge_aqueous_pressure = merge_aq; if (transport_conf.has_attribute(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_CUTOFFRESIDUAL)) { opts.cutoff_residuals = transport_conf.get_attribute( SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_CUTOFFRESIDUAL); } SPC_CONF_LOG << "Cutoff |R|_0^2 : " << opts.cutoff_residuals; auto stagger = UnsaturatedTransportStagger::make( simul_data.variables, simul_data.bcs, opts); // ! 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.str(std::string()); 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.str(std::string()); 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_gas_options(aq)) = opts; } } SPC_CONF_LOG << "\nTransport options for gaseous equations \n" << SPC_CONF_LOG_HLINE; msg.str(std::string()); 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_gas_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); if (not simul_data.variables->component_has_gas(comp_id)) { std::string msg = "Component " + comp_label + " has no gas." \ "Can't parse the configuration for gas equation solver"; ERROR << msg; throw std::runtime_error(msg); } // H2O if (comp_id == 0) { if (merge_sat) continue; io::configure_transport_options( *(stagger->get_gas_options(0)), std::move(subconf) ); SPC_CONF_LOG << "\nTransport options for gaseous equations " << comp_label << ". \n" << SPC_CONF_LOG_HLINE; msg.str(std::string()); io::print_transport_options(*(stagger->get_gas_options(0)), msg); SPC_CONF_LOG << msg.str() << SPC_CONF_LOG_HLINE; } // aqueous component else { if (merge_aq) continue; io::configure_transport_options( *(stagger->get_gas_options(comp_id)), std::move(subconf) ); SPC_CONF_LOG << "\nTransport options for gaseous equations " << comp_label << ". \n" << SPC_CONF_LOG_HLINE; msg.str(std::string()); io::print_transport_options(*(stagger->get_gas_options(comp_id)), msg); SPC_CONF_LOG << msg.str() << SPC_CONF_LOG_HLINE; } } } } return stagger; } std::string initialize_output( SimulationData& simul_data, const std::string& work_dir, io::YAMLConfigHandle&& conf, io::AllIOFiles& all_io_file ) { 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); } all_io_file.add_solution(io::output_file(IO_FILE_OUTPUT, 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); all_io_file.add_database(io::output_file(IO_FILE_WORK_DB, 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; }