diff --git a/examples/reactmicp/saturated_react/momas_benchmark.cpp b/examples/reactmicp/saturated_react/momas_benchmark.cpp index 2e67c53..73fbe3c 100644 --- a/examples/reactmicp/saturated_react/momas_benchmark.cpp +++ b/examples/reactmicp/saturated_react/momas_benchmark.cpp @@ -1,594 +1,594 @@ /* ============================================================================= 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 "reactmicp/reactmicp.hpp" #include "specmicp_database/aqueous_selector.hpp" #include #include #include #include // parameters #define DX 0.05 #define IS_ADVECTIVE true #define MIN_DT 0.001 #define MAX_DT 100.0 #define RESTART_DT MIN_DT using namespace specmicp; RawDatabasePtr get_momas_db() { database::Database thedatabase("../data/momas_benchmark.yaml"); return thedatabase.get_database(); } void prepare_db_for_easy_case(RawDatabasePtr raw_data) { database::Database db_handler(raw_data); db_handler.remove_components({"X5",}); if (raw_data->nb_component() != 6) { throw std::runtime_error("Not enough or Too many component"); } db_handler.remove_solid_phases(); if (not raw_data->is_valid()) { throw std::runtime_error("Database is not valid."); } //std::cout << selector.aqueous_label_to_id("C6") << " - " << selector.aqueous_label_to_id("C7") << std::endl; database::AqueousSelector selector(raw_data); selector.remove_aqueous({selector.aqueous_label_to_id("C6"), selector.aqueous_label_to_id("C7")}); std::cout << raw_data->aqueous.get_nu_matrix() << std::endl; if (raw_data->nb_aqueous() != 5) { throw std::runtime_error("Not enough or Too many aqueous."); } if (not raw_data->is_valid()) { throw std::runtime_error("Database is not valid."); } } AdimensionalSystemSolution easy_material_a( RawDatabasePtr raw_data, const AdimensionalSystemConstraints& constraints, // const AdimensionalSystemSolution& initial_solution, const AdimensionalSystemSolverOptions& options ) { specmicp::Vector variables ; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, true); if (current_perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized) throw std::runtime_error("Failed to solve initial solution for material A."); return solver.get_raw_solution(variables); } AdimensionalSystemSolution easy_material_b( RawDatabasePtr raw_data, const AdimensionalSystemConstraints& constraints, const specmicp::AdimensionalSystemSolverOptions& options ) { specmicp::Vector variables; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, true); if (current_perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized) throw std::runtime_error("Failed to solve initial solution for material B."); return solver.get_raw_solution(variables); } AdimensionalSystemSolution easy_bc( RawDatabasePtr raw_data, const AdimensionalSystemConstraints& constraints, const specmicp::AdimensionalSystemSolverOptions& options ) { specmicp::Vector variables; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); - solver.initialise_variables(variables, 1.0, -4.0); + solver.initialize_variables(variables, 1.0, -4.0); solver.run_pcfm(variables); specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, false); if (current_perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized) throw std::runtime_error("Failed to solve initial solution for BC."); return solver.get_raw_solution(variables); } using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::solver; using namespace specmicp::reactmicp::systems::satdiff; class MoMasUspcaling: public solver::UpscalingStaggerBase { public: MoMasUspcaling( RawDatabasePtr the_database, units::UnitsSet the_units, bool advective ): m_data(the_database), m_units_set(the_units), m_advective(advective) { } //! \brief Initialize the stagger at the beginning of the computation void initialize(VariablesBasePtr var) { return initialize(var.get()); } virtual void initialize(VariablesBase * const var) override { SaturatedVariables * const true_var = static_cast(var); database::Database db_handler(m_data); true_var->upscaling_variables().setZero(); for (index_t node=0; nodenb_nodes(); ++node) { // if material b const scalar_t coord = true_var->get_mesh()->get_position(node); if (coord < 1.1 and coord >= 1.0) { true_var->porosity(node) = 0.5; true_var->permeability(node) = 1e-5; if (m_advective) true_var->diffusion_coefficient(node) = 3.3e-4; else true_var->diffusion_coefficient(node) = 330e-3; true_var->fluid_velocity(node) = 5.5e-3; } else // material a { true_var->porosity(node) = 0.25; true_var->permeability(node) = 1e-2; if (m_advective) true_var->diffusion_coefficient(node) = 5.5e-5; else true_var->diffusion_coefficient(node) = 55e-3; true_var->fluid_velocity(node) = 5.5e-3; } std::cout << "node : " << node << " - porosity : " << true_var->porosity(node) << std::endl; } } //! \brief Initialize the stagger at the beginning of an iteration virtual void initialize_timestep(scalar_t _, VariablesBase * const __) override { } //! \brief Solve the equation for the timestep virtual StaggerReturnCode restart_timestep(VariablesBase * const _) override { return StaggerReturnCode::ResidualMinimized; } private: RawDatabasePtr m_data; units::UnitsSet m_units_set; bool m_advective; }; AdimensionalSystemConstraints get_constraints_easy_a(RawDatabasePtr& raw_data) { Vector total_concentrations(Vector::Zero(raw_data->nb_component())); total_concentrations(raw_data->get_id_component("X2")) = -2.0; total_concentrations(raw_data->get_id_component("X4")) = 2.0; AdimensionalSystemConstraints constraints_a(0.25*total_concentrations); constraints_a.disable_conservation_water(); constraints_a.surface_model.model_type = SurfaceEquationType::Equilibrium; constraints_a.surface_model.concentration = 0.25*1.0; constraints_a.inert_volume_fraction = 0.75; return constraints_a; } AdimensionalSystemConstraints get_constraints_easy_b(RawDatabasePtr& raw_data) { Vector total_concentrations(Vector::Zero(raw_data->nb_component())); total_concentrations(raw_data->get_id_component("X2")) = -2.0; total_concentrations(raw_data->get_id_component("X4")) = 2.0; AdimensionalSystemConstraints constraints_b(0.5*total_concentrations); constraints_b.disable_conservation_water(); constraints_b.surface_model.model_type = SurfaceEquationType::Equilibrium; constraints_b.surface_model.concentration = 0.5*10.0; constraints_b.inert_volume_fraction = 0.5; return constraints_b; } AdimensionalSystemConstraints get_constraints_easy_bc(RawDatabasePtr& raw_data) { Vector total_concentrations(Vector::Zero(raw_data->nb_component())); total_concentrations(raw_data->get_id_component("X1")) = 0.3; total_concentrations(raw_data->get_id_component("X2")) = 0.3; total_concentrations(raw_data->get_id_component("X3")) = 0.3; AdimensionalSystemConstraints constraints_bc(total_concentrations); constraints_bc.disable_conservation_water(); constraints_bc.surface_model.model_type = SurfaceEquationType::NoEquation; constraints_bc.surface_model.concentration = 0.0; constraints_bc.inert_volume_fraction = 0.0; return constraints_bc; } AdimensionalSystemSolverOptions get_specmicp_options() { AdimensionalSystemSolverOptions options; options.solver_options.maxstep = 1.0; options.solver_options.max_iter = 200; options.solver_options.maxiter_maxstep = 200; options.solver_options.use_crashing = false; options.solver_options.use_scaling = false; options.solver_options.non_monotone_linesearch = false; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 100; options.solver_options.projection_min_variable = 1e-9; options.solver_options.fvectol = 1e-10; options.solver_options.steptol = 1e-14; options.solver_options.threshold_cycling_linesearch = 1e-8; options.system_options.non_ideality_tolerance = 1e-10; options.system_options.non_ideality = false; options.system_options.restart_concentration = -2; options.use_pcfm = true; return options; } mesh::Mesh1DPtr get_mesh(scalar_t dx) { scalar_t tot_length = 2.1+dx; index_t nb_nodes = tot_length/dx +2; return mesh::uniform_mesh1d(nb_nodes, dx, 1.0); } void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options) { transport_options.maximum_iterations = 450; transport_options.quasi_newton = 3; transport_options.residuals_tolerance = 1e-8; transport_options.step_tolerance = 1e-12; transport_options.absolute_tolerance = 1e-6; transport_options.alpha = 1.0; //transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang; transport_options.sparse_solver = SparseSolver::SparseLU; //transport_options.sparse_solver = SparseSolver::GMRES; transport_options.threshold_stationary_point = 1e-2; } void set_reactive_solver_options(reactmicp::solver::ReactiveTransportOptions& reactive_options) { reactive_options.maximum_iterations = 50; //500; reactive_options.residuals_tolerance = 1e-2; reactive_options.good_enough_tolerance = 1.0; reactive_options.absolute_residuals_tolerance = 1e-6; reactive_options.implicit_upscaling = false; } int main() { Eigen::initParallel(); Timer global_timer; global_timer.start(); io::OutFile output_gen("out_momas_out.txt"); io::print_reactmicp_header(output_gen); specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; RawDatabasePtr database = get_momas_db(); prepare_db_for_easy_case(database); units::UnitsSet the_units = units::UnitsSet(); the_units.length = units::LengthUnit::decimeter; // so that rho_w ~ 1 scalar_t dx=DX; mesh::Mesh1DPtr the_mesh = get_mesh(dx); AdimensionalSystemConstraints constraints_a = get_constraints_easy_a(database); AdimensionalSystemConstraints constraints_b = get_constraints_easy_b(database); AdimensionalSystemConstraints constraints_bc = get_constraints_easy_bc(database); AdimensionalSystemSolverOptions options = get_specmicp_options(); options.units_set = the_units; AdimensionalSystemSolution mat_bc = easy_bc(database, constraints_bc, options); AdimensionalSystemSolution mat_a = easy_material_a(database, constraints_a, options); AdimensionalSystemSolution mat_b = easy_material_b(database, constraints_b, options); options.solver_options.maxstep = 1.0; options.solver_options.use_scaling = true; options.solver_options.non_monotone_linesearch = true; index_t nb_nodes = the_mesh->nb_nodes(); std::vector list_constraints = {constraints_a, constraints_b, constraints_bc}; std::vector list_initial_composition = {mat_a, mat_b, mat_bc}; std::vector index_initial_state(nb_nodes, 0); std::vector index_constraints(nb_nodes, 0); index_t start_mat_b = 1.0/dx; index_t end_mat_b = 1.1/dx; for (index_t node=start_mat_b; node list_fixed_nodes = {0, }; std::map list_slaves_nodes = {{nb_nodes-1, nb_nodes-2}}; std::vector list_immobile_components = {0, 1}; systems::satdiff::SaturatedVariablesPtr variables = systems::satdiff::init_variables(the_mesh, database, the_units, list_fixed_nodes, list_initial_composition, index_initial_state); ChemistryStaggerPtr equilibrium_stagger = std::make_shared(list_constraints, index_constraints, options); std::shared_ptr transport_stagger = std::make_shared( variables, list_fixed_nodes, list_slaves_nodes, list_immobile_components); set_transport_options(transport_stagger->options_solver()); const bool is_advective = IS_ADVECTIVE; UpscalingStaggerPtr upscaling_stagger = std::make_shared(database, the_units, is_advective); ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger); transport_stagger->initialize(variables); equilibrium_stagger->initialize(variables); upscaling_stagger->initialize(variables); set_reactive_solver_options(solver.get_options()); // mesh output io::print_mesh("out_momas_mesh.dat", the_mesh, the_units.length); //std::cout << variables->displacement() << std::endl; io::OutFile output_iter("out_momas_iter.dat"); std::ofstream output_x1, output_x2, output_x3, output_x4; std::ofstream output_sx1, output_sx2, output_sx3, output_sx4; std::ofstream output_c1, output_c2, output_s; output_x1.open("out_momas_x1.dat"); output_x2.open("out_momas_x2.dat"); output_x3.open("out_momas_x3.dat"); output_x4.open("out_momas_x4.dat"); output_sx1.open("out_momas_sx1.dat"); output_sx2.open("out_momas_sx2.dat"); output_sx3.open("out_momas_sx3.dat"); output_sx4.open("out_momas_sx4.dat"); output_c1.open("out_momas_c1.dat"); output_c2.open("out_momas_c2.dat"); output_s.open("out_momas_s.dat"); index_t cnt = 0; scalar_t dt=RESTART_DT; reactmicp::solver::Timestepper timestepper(MIN_DT, MAX_DT, 5000, 2.0); timestepper.get_options().alpha_average = 0.3; timestepper.get_options().decrease_failure = 0.1; timestepper.get_options().increase_factor = 1.05; timestepper.get_options().decrease_factor = 0.50; timestepper.get_options().iteration_lower_target = 1.001; io::print_reactmicp_performance_long_header(output_iter); while (timestepper.get_total() < timestepper.get_total_target()) { output_gen << "dt : " << dt << std::endl; Timer step_timer; step_timer.start(); reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables); step_timer.stop(); // if (retcode > reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) //{ io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total()+dt, solver.get_perfs()); //} dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations); if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) { dt = RESTART_DT; variables->reset_main_variables(); retcode = solver.solve_timestep(dt, variables); dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations); io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total(), solver.get_perfs()); } output_gen << "Time : " << timestepper.get_total() << std::endl; io::print_reactmicp_performance_short(output_gen, solver.get_perfs()); scalar_t total = timestepper.get_total(); if (cnt % 10== 0) { output_x1 << total; output_x2 << total; output_x3 << total; output_x4 << total; output_sx1 << total; output_sx2 << total; output_sx3 << total; output_sx4 << total; output_c1 << total; output_c2 << total; output_s << total; for (index_t node: the_mesh->range_nodes()) { output_x1 << "\t" << variables->aqueous_concentration(node, 1, variables->displacement()); output_x2 << "\t" << variables->aqueous_concentration(node, 2, variables->displacement()); output_x3 << "\t" << variables->aqueous_concentration(node, 3, variables->displacement()); output_x4 << "\t" << variables->aqueous_concentration(node, 4, variables->displacement()); output_sx1 << "\t" << variables->solid_concentration(node, 1, variables->displacement()); output_sx2 << "\t" << variables->solid_concentration(node, 2, variables->displacement()); output_sx3 << "\t" << variables->solid_concentration(node, 3, variables->displacement()); output_sx4 << "\t" << variables->solid_concentration(node, 4, variables->displacement()); output_c1 << "\t" << variables->equilibrium_solution(node).secondary_molalities(0); output_c2 << "\t" << variables->equilibrium_solution(node).secondary_molalities(1); output_s << "\t" << pow10(variables->equilibrium_solution(node).main_variables(5)); } output_x1 << std::endl; output_x2 << std::endl; output_x3 << std::endl; output_x4 << std::endl; output_sx1 << std::endl; output_sx2 << std::endl; output_sx3 << std::endl; output_sx4 << std::endl; output_c1 << std::endl; output_c2 << std::endl; output_s << std::endl; } ++ cnt; } variables->displacement()(database->get_id_component("X1")) = 0; variables->displacement()(database->get_id_component("X2")) = -2.0; variables->displacement()(database->get_id_component("X3")) = 0; variables->displacement()(database->get_id_component("X4")) = 2.0; timestepper.set_total_target(6000); dt = RESTART_DT; transport_stagger->options_solver().absolute_tolerance = 1e-10; solver.get_options().absolute_residuals_tolerance = 1e-6; while (timestepper.get_total() < timestepper.get_total_target()) { output_gen << "dt : " << dt << std::endl; Timer step_timer; step_timer.start(); reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables); step_timer.stop(); // if (retcode > reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) //{ io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total()+dt, solver.get_perfs()); //} dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations); if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) { dt = RESTART_DT; variables->reset_main_variables(); retcode = solver.solve_timestep(dt, variables); dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations); io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total(), solver.get_perfs()); } output_gen << "Time : " << timestepper.get_total() << std::endl; io::print_reactmicp_performance_short(output_gen, solver.get_perfs()); scalar_t total = timestepper.get_total(); if (cnt % 10== 0) { output_x1 << total; output_x2 << total; output_x3 << total; output_x4 << total; output_sx1 << total; output_sx2 << total; output_sx3 << total; output_sx4 << total; output_c1 << total; output_c2 << total; output_s << total; for (index_t node: the_mesh->range_nodes()) { output_x1 << "\t" << variables->aqueous_concentration(node, 1, variables->displacement()); output_x2 << "\t" << variables->aqueous_concentration(node, 2, variables->displacement()); output_x3 << "\t" << variables->aqueous_concentration(node, 3, variables->displacement()); output_x4 << "\t" << variables->aqueous_concentration(node, 4, variables->displacement()); output_sx1 << "\t" << variables->solid_concentration(node, 1, variables->displacement()); output_sx2 << "\t" << variables->solid_concentration(node, 2, variables->displacement()); output_sx3 << "\t" << variables->solid_concentration(node, 3, variables->displacement()); output_sx4 << "\t" << variables->solid_concentration(node, 4, variables->displacement()); output_c1 << "\t" << variables->equilibrium_solution(node).secondary_molalities(0); output_c2 << "\t" << variables->equilibrium_solution(node).secondary_molalities(1); output_s << "\t" << pow10(variables->equilibrium_solution(node).main_variables(5)); } output_x1 << std::endl; output_x2 << std::endl; output_x3 << std::endl; output_x4 << std::endl; output_sx1 << std::endl; output_sx2 << std::endl; output_sx3 << std::endl; output_sx4 << std::endl; output_c1 << std::endl; output_c2 << std::endl; output_s << std::endl; } ++ cnt; } reactmicp::solver::ReactiveTransportTimer timer = solver.get_timer(); io::print_reactmicp_timer(output_gen, timer); global_timer.stop(); output_gen << "Total computation time : " << global_timer.elapsed_time() << "s." << std::endl; return EXIT_SUCCESS; } diff --git a/examples/reactmicp/saturated_react/react_leaching.cpp b/examples/reactmicp/saturated_react/react_leaching.cpp index b32c1c9..4eb4205 100644 --- a/examples/reactmicp/saturated_react/react_leaching.cpp +++ b/examples/reactmicp/saturated_react/react_leaching.cpp @@ -1,764 +1,764 @@ /* ============================================================================= 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 "reactmicp/reactmicp.hpp" #include #include #include using namespace specmicp; // Initial conditions // ================== specmicp::RawDatabasePtr leaching_database() { specmicp::database::Database thedatabase("../data/cemdata_specmicp.yaml"); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"}, {"Al[3+]","Al(OH)4[-]"} }); thedatabase.remove_gas_phases(); thedatabase.swap_components(swapping); return thedatabase.get_database(); } const double M_CaO = 56.08; const double M_SiO2 = 60.09; const double M_Al2O3 = 101.96; const double M_SO3 = 80.06; void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species) { 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::FullPivLU solver(compo_in_oxyde); compo_species = solver.solve(compo_oxyde); } AdimensionalSystemSolution initial_leaching_pb( RawDatabasePtr raw_data, scalar_t mult, units::UnitsSet the_units) { Formulation formulation; Vector oxyde_compo(4); oxyde_compo << 62.9, 20.6, 5.8, 3.1; Vector species_compo; compo_from_oxyde(oxyde_compo, species_compo); std::cout << species_compo << std::endl; scalar_t m_c3s = mult*species_compo(0); //0.6; scalar_t m_c2s = mult*species_compo(1); //0.2; scalar_t m_c3a = mult*species_compo(2); //0.1; scalar_t m_gypsum = mult*species_compo(3); //0.1; scalar_t wc = 0.5; scalar_t m_water = wc*1.0e-3*( m_c3s*(3*56.08+60.08) + m_c2s*(2*56.06+60.08) + m_c3a*(3*56.08+101.96) + m_gypsum*(56.08+80.06+2*18.02) ); formulation.mass_solution = m_water; formulation.amount_minerals = { {"C3S", m_c3s}, {"C2S", m_c2s}, {"C3A", m_c3a}, {"Gypsum", m_gypsum} }; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = 1; constraints.set_saturated_system(); specmicp::AdimensionalSystemSolverOptions options; options.solver_options.maxstep = 100.0; options.solver_options.max_iter = 200; options.solver_options.maxiter_maxstep = 100; options.solver_options.projection_min_variable = 1e-9; options.solver_options.set_tolerance(1e-10, 1e-14); options.solver_options.disable_crashing(); options.solver_options.enable_scaling(); options.solver_options.disable_descent_direction(); options.solver_options.enable_non_monotone_linesearch(); options.system_options.non_ideality_tolerance = 1e-10; options.units_set = the_units; specmicp::Vector x; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); - solver.initialise_variables(x, 0.8, -3.0); + solver.initialize_variables(x, 0.8, -3.0); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false); std::cout << "Initial return code : " << (int) perf.return_code << std::endl; if ((int) perf.return_code <= (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet) { throw std::runtime_error("Initial solution has not converged"); } return solver.get_raw_solution(x); } AdimensionalSystemSolution initial_blank_leaching_pb( RawDatabasePtr raw_data, scalar_t mult, units::UnitsSet the_units) { Formulation formulation; Vector oxyde_compo(4); oxyde_compo << 62.9, 20.6, 5.8, 3.1; Vector species_compo; compo_from_oxyde(oxyde_compo, species_compo); std::cout << species_compo << std::endl; scalar_t m_c3s = 0.6; scalar_t m_c2s = 0.2; scalar_t m_c3a = 0.1; scalar_t m_gypsum = 0.1; scalar_t wc = 0.8; scalar_t m_water = wc*1.0e-3*( m_c3s*(3*56.08+60.08) + m_c2s*(2*56.06+60.08) + m_c3a*(3*56.08+101.96) + m_gypsum*(56.08+80.06+2*18.02) ); formulation.mass_solution = m_water; formulation.amount_minerals = { {"SiO2,am", mult}, }; // formulation.concentration_aqueous = { // {"Ca[2+]", 1e-7}, // {"Al(OH)4[-]", 1e-7}, // {"SO4[2-]", 1e-7} // }; formulation.extra_components_to_keep = {"Ca[2+]", "SO4[2-]", "Al(OH)4[-]"}; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = 1; constraints.set_saturated_system(); specmicp::AdimensionalSystemSolverOptions options; options.solver_options.maxstep = 10.0; options.solver_options.max_iter = 200; options.solver_options.maxiter_maxstep = 200; options.solver_options.use_crashing = false; options.solver_options.use_scaling = true; options.solver_options.factor_descent_condition = -1.0; options.solver_options.factor_gradient_search_direction = 100; options.solver_options.projection_min_variable = 1e-9; options.solver_options.fvectol = 1e-4; options.solver_options.steptol = 1e-14; options.solver_options.non_monotone_linesearch = true; options.system_options.non_ideality_tolerance = 1e-10; options.units_set = the_units; Vector x; AdimensionalSystemSolver solver(raw_data, constraints, options); - solver.initialise_variables(x, 0.75, -4.0); + solver.initialize_variables(x, 0.75, -4.0); micpsolver::MiCPPerformance perf = solver.solve(x, false); std::cout << "Initial return code : " << (int) perf.return_code << std::endl; if ((int) perf.return_code <= (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet) { throw std::runtime_error("Initial solution has not converged"); } AdimensionalSystemSolution solution = solver.get_raw_solution(x); AdimensionalSystemSolutionModificator extractor(solution, raw_data, the_units); extractor.remove_solids(); Vector totconc(raw_data->nb_component()); totconc(0) = extractor.density_water()*extractor.total_aqueous_concentration(0); for (index_t component: raw_data->range_aqueous_component()) { totconc(component) = 0.01*extractor.density_water()*extractor.total_aqueous_concentration(component); } constraints.total_concentrations = totconc; solver = AdimensionalSystemSolver(raw_data, constraints, options); perf = solver.solve(x, false); specmicp_assert((int) perf.return_code > (int) micpsolver::MiCPSolverReturnCode::NotConvergedYet); return solver.get_raw_solution(x); } // Upscaling // ========= using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::solver; using namespace specmicp::reactmicp::systems::satdiff; class LeachingUspcaling: public solver::UpscalingStaggerBase { public: LeachingUspcaling(RawDatabasePtr the_database, units::UnitsSet the_units): m_data(the_database), m_units_set(the_units) { } //! \brief Initialize the stagger at the beginning of the computation void initialize(VariablesBasePtr var) { return initialize(var.get()); } virtual void initialize(VariablesBase * const var) override { SaturatedVariables * const true_var = static_cast(var); database::Database db_handler(m_data); m_dt = 1; m_id_cshj = db_handler.mineral_label_to_id("CSH,j"); m_initial_sat_cshj = true_var->equilibrium_solution(1).main_variables( m_data->nb_component()+m_id_cshj); //m_initial_sat_cshj = 0.36; std::cout << m_initial_sat_cshj << std::endl; true_var->upscaling_variables().setZero(); for (index_t node=0; nodenb_nodes(); ++node) { upscaling_one_node(node, true_var); } } //! \brief Initialize the stagger at the beginning of an iteration virtual void initialize_timestep(scalar_t dt, VariablesBase * const var) override { SaturatedVariables * const true_var = static_cast(var); m_dt = dt; for (index_t node=1; nodenb_nodes(); ++node) { upscaling_one_node(node, true_var); true_var->vel_porosity(node) = 0.0; } } //! \brief Solve the equation for the timestep virtual StaggerReturnCode restart_timestep(VariablesBase * const var) override { SaturatedVariables* true_var = static_cast(var); for (index_t node=1; nodenb_nodes(); ++node) { upscaling_one_node(node, true_var); } return StaggerReturnCode::ResidualMinimized; } scalar_t diffusion_coefficient(scalar_t porosity) { const scalar_t factor = 1e4; scalar_t poro = factor*std::exp(9.85*porosity-29.08); // const scalar_t de_0 = 2.76e-12; // const scalar_t porosity_r = 0.02; // const scalar_t exponent = 3.32; //if (porosity > 0.7) // return factor*4.0*1e-10; // else // { // const scalar_t ratio = ((porosity - porosity_r)/(0.25 - porosity_r)); // poro = factor*de_0*std::pow(ratio, exponent); // } return std::min(poro,factor*1e-10); } void upscaling_one_node(index_t node, SaturatedVariables * const true_var) { AdimensionalSystemSolutionExtractor extractor( true_var->equilibrium_solution(node), m_data, m_units_set ); scalar_t porosity = extractor.porosity(); //- true_var->equilibrium_solution(node).main_variables(m_data->nb_component)/m_initial_sat_cshj*0.05; //std::cout << "porosity : " << porosity // << " - saturation water" << true_var->equilibrium_solution(node).main_variables(0) << std::endl; true_var->vel_porosity(node) += (porosity - true_var->porosity(node))/m_dt; true_var->porosity(node) = porosity; true_var->diffusion_coefficient(node) = diffusion_coefficient(porosity); } private: RawDatabasePtr m_data; units::UnitsSet m_units_set; index_t m_id_cshj; scalar_t m_initial_sat_cshj; scalar_t m_dt; }; // void set_specmicp_options(AdimensionalSystemSolverOptions& options) { options.solver_options.maxstep = 10.0; options.solver_options.max_iter = 100; options.solver_options.maxiter_maxstep = 100; options.solver_options.use_crashing = false; options.solver_options.use_scaling = true; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 100; options.solver_options.projection_min_variable = 1e-16; options.solver_options.fvectol = 1e-10; options.solver_options.steptol = 1e-14; options.solver_options.non_monotone_linesearch = true; options.system_options.non_ideality_tolerance = 1e-12; options.system_options.non_ideality_max_iter = 100; options.solver_options.threshold_cycling_linesearch = 1e-5; } void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options) { transport_options.maximum_iterations = 450; transport_options.quasi_newton = 3; transport_options.residuals_tolerance = 1e-6; transport_options.step_tolerance = 1e-14; transport_options.alpha = 1; transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang; transport_options.sparse_solver = SparseSolver::SparseLU; //transport_options.sparse_solver = SparseSolver::GMRES; transport_options.threshold_stationary_point = 1e-2; transport_options.absolute_tolerance = 1e-16; } void set_reactive_solver_options(ReactiveTransportOptions& solver_options) { solver_options.maximum_iterations = 100; solver_options.residuals_tolerance = 1e-2; solver_options.good_enough_tolerance = 1.0; solver_options.absolute_residuals_tolerance = 1e-16; solver_options.implicit_upscaling = true; } mesh::Mesh1DPtr get_mesh(scalar_t dx, index_t nb_nodes, index_t additional_nodes) { mesh::UniformAxisymmetricMesh1DGeometry geometry; geometry.dx = dx; geometry.radius = 3.5 + additional_nodes*dx; //cm geometry.height = 8.0; geometry.nb_nodes = nb_nodes + additional_nodes; return mesh::uniform_axisymmetric_mesh1d(geometry); } mesh::Mesh1DPtr get_mesh() { Vector radius(66); const scalar_t height = 8.0; radius << 3.5005, 3.5000, 3.4995, 3.4990, 3.4980, 3.4970, 3.4960, 3.4950, 3.4925, 3.4900, 3.4875, 3.4850, 3.4825, 3.4800, 3.4775, 3.4750, 3.4725, 3.4700, 3.4675, 3.4650, 3.4625, 3.4600, 3.4575, 3.4550, 3.4525, 3.4500, 3.4475, 3.445, 3.4425, 3.440, 3.4375, 3.435, 3.4325, 3.430, 3.4275, 3.425, 3.4225, 3.420, 3.4175, 3.415, 3.4125, 3.410, 3.4075, 3.405, 3.4025, 3.400, 3.3975, 3.395, 3.3925, 3.390, 3.3875, 3.385, 3.3825, 3.38, 3.3775, 3.3750, 3.3725, 3.3700, 3.3675, 3.365, 3.3625, 3.36, 3.3575, 3.355, 3.325, 3.3 ; return mesh::axisymmetric_mesh1d(radius, height); } void output_mesh(mesh::Mesh1DPtr the_mesh) { std::ofstream output_mesh; output_mesh.open("out_leaching_mesh.dat"); output_mesh << "Node\tPosition" << std::endl; for (index_t node: the_mesh->range_nodes()) { output_mesh << node << "\t" << the_mesh->get_position(node) << std::endl; } output_mesh.close(); } void print_minerals_profile( RawDatabasePtr the_database, SaturatedVariablesPtr variables, mesh::Mesh1DPtr the_mesh, const std::string& filepath ) { std::ofstream ofile(filepath); ofile << "Radius"; for (index_t mineral: the_database->range_mineral()) { ofile << "\t" << the_database->get_label_mineral(mineral); } ofile << "\t" << "Porosity"; ofile << "\t" << "pH"; ofile << std::endl; for (index_t node: the_mesh->range_nodes()) { ofile << the_mesh->get_position(node); AdimensionalSystemSolutionExtractor extractor(variables->equilibrium_solution(node), the_database, units::UnitsSet()); for (index_t mineral: the_database->range_mineral()) { ofile << "\t" << extractor.volume_fraction_mineral(mineral); } ofile << "\t" << variables->porosity(node); ofile << "\t" << extractor.pH(); ofile << std::endl; } ofile.close(); } // Main // ==== int main() { Eigen::initParallel(); specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; RawDatabasePtr raw_data = leaching_database(); units::UnitsSet the_units; the_units.length = units::LengthUnit::centimeter; AdimensionalSystemSolution initial = initial_leaching_pb(raw_data, 0.015, the_units); AdimensionalSystemSolutionModificator extractor(initial, raw_data, the_units); index_t id_ca = raw_data->get_id_component("Ca[2+]"); index_t id_si = raw_data->get_id_component("SiO(OH)3[-]"); std::cout << extractor.total_solid_concentration(id_ca) << " - " << extractor.porosity() << std::endl; extractor.scale_total_concentration(id_ca, 0.0145); std::cout << extractor.total_solid_concentration(id_ca) << " - " << extractor.porosity() << std::endl; std::cout << initial.main_variables << std::endl; AdimensionalSystemSolution blank = initial_blank_leaching_pb(raw_data, 0.005, the_units); std::cout << blank.main_variables << std::endl; scalar_t dx = 0.0005; index_t nb_nodes = 150; index_t additional_nodes = 1; //mesh::Mesh1DPtr the_mesh = get_mesh(dx, nb_nodes, additional_nodes); mesh::Mesh1DPtr the_mesh = get_mesh(); nb_nodes = the_mesh->nb_nodes(); std::vector list_initial_composition = {initial, blank}; std::vector index_initial_state(nb_nodes, 0); for (int i=0; i< additional_nodes; ++i) index_initial_state[i] = 1; std::vector list_fixed_nodes = {0, }; systems::satdiff::SaturatedVariablesPtr variables = systems::satdiff::init_variables(the_mesh, raw_data, the_units, list_fixed_nodes, list_initial_composition, index_initial_state); specmicp::AdimensionalSystemConstraints constraints; constraints.charge_keeper = 1; constraints.set_saturated_system(); constraints.inert_volume_fraction = 0.0; AdimensionalSystemSolverOptions options; set_specmicp_options(options); options.units_set = the_units; ChemistryStaggerPtr equilibrium_stagger = std::make_shared(nb_nodes, constraints, options); std::shared_ptr transport_stagger = std::make_shared(variables, list_fixed_nodes); dfpmsolver::ParabolicDriverOptions& transport_options = transport_stagger->options_solver(); set_transport_options(transport_options); UpscalingStaggerPtr upscaling_stagger = std::make_shared(raw_data, the_units); ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger); transport_stagger->initialize(variables); equilibrium_stagger->initialize(variables); upscaling_stagger->initialize(variables); set_reactive_solver_options(solver.get_options()); io::OutFile output_iter("out_leaching_iter.dat"); std::ofstream output_c_ca, output_s_ca, output_s_si; std::ofstream output_poro; std::ofstream output_loss_ca; output_c_ca.open("out_c_ca.dat"); output_s_ca.open("out_s_ca.dat"); output_s_si.open("out_s_si.dat"); output_poro.open("out_poro.dat"); output_loss_ca.open("out_loss_ca.dat"); scalar_t initial_ca = 0; for (index_t node: the_mesh->range_nodes()) { initial_ca += variables->solid_concentration(node, id_ca, variables->displacement()) * the_mesh->get_volume_cell(node); } scalar_t dt=2.0; io::print_reactmicp_performance_long_header(output_iter); Timestepper timestepper(1.0, 500.0, 25*24*3600, 2); int cnt =0; while (timestepper.get_total() < timestepper.get_total_target()) { Timer step_timer; step_timer.start(); reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables); step_timer.stop(); io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total() , solver.get_perfs()); dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations); if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) { dt = 1.0; variables->reset_main_variables(); retcode = solver.solve_timestep(dt, variables); dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations); io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total(), solver.get_perfs()); } // std::cout << "time : " << total/3600/24 << std::endl // << "Iter : " << solver.get_perfs().nb_iterations << std::endl // << "Residuals : " << solver.get_perfs().residuals << std::endl; if (cnt % 100 == 0) { scalar_t time = timestepper.get_total()/3600.0/24.0; scalar_t sqrt_time = std::sqrt(time); output_c_ca << time; output_s_ca << time; output_s_si << time; output_loss_ca << time << "\t" << sqrt_time; output_poro << time; scalar_t amount_ca = 0.0; for (index_t node: the_mesh->range_nodes()) { amount_ca += variables->solid_concentration(node, id_ca, variables->displacement()) * the_mesh->get_volume_cell(node); output_c_ca << "\t" << variables->aqueous_concentration(node, id_ca, variables->displacement()); output_s_ca << "\t" << variables->solid_concentration( node, id_ca, variables->displacement()); output_s_si << "\t" << variables->solid_concentration( node, id_si, variables->displacement()); output_poro << "\t" << variables->porosity(node); } output_poro << std::endl; output_loss_ca << "\t" << (initial_ca - amount_ca)/1.75929e-2 << std::endl; output_c_ca << std::endl; output_s_ca << std::endl; output_s_si << std::endl; } ++cnt; } print_minerals_profile(raw_data, variables, the_mesh, "out_leaching_minerals_profile_25d.dat"); timestepper.set_total_target(90*24*3600); while (timestepper.get_total() < timestepper.get_total_target()) { Timer step_timer; step_timer.start(); reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables); step_timer.stop(); io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total() , solver.get_perfs()); dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations); if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet) { dt = 1.0; variables->reset_main_variables(); retcode = solver.solve_timestep(dt, variables); dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations); io::print_reactmicp_performance_long(output_iter, cnt, timestepper.get_total(), solver.get_perfs()); } // std::cout << "time : " << total/3600/24 << std::endl // << "Iter : " << solver.get_perfs().nb_iterations << std::endl // << "Residuals : " << solver.get_perfs().residuals << std::endl; if (cnt % 100 == 0) { scalar_t time = timestepper.get_total()/3600.0/24.0; scalar_t sqrt_time = std::sqrt(time); output_c_ca << time; output_s_ca << time; output_s_si << time; output_loss_ca << time << "\t" << sqrt_time; output_poro << time; scalar_t amount_ca = 0.0; for (index_t node: the_mesh->range_nodes()) { amount_ca += variables->solid_concentration(node, id_ca, variables->displacement()) * the_mesh->get_volume_cell(node); output_c_ca << "\t" << variables->aqueous_concentration(node, id_ca, variables->displacement()); output_s_ca << "\t" << variables->solid_concentration( node, id_ca, variables->displacement()); output_s_si << "\t" << variables->solid_concentration( node, id_si, variables->displacement()); output_poro << "\t" << variables->porosity(node); } output_poro << std::endl; output_loss_ca << "\t" << (initial_ca - amount_ca)/1.75929e-2 << std::endl; output_c_ca << std::endl; output_s_ca << std::endl; output_s_si << std::endl; } ++cnt; } io::print_reactmicp_timer(output_iter, solver.get_timer()); } diff --git a/examples/reactmicp/unsaturated/acc_carbo.cpp b/examples/reactmicp/unsaturated/acc_carbo.cpp index 5b95d45..f7480b3 100644 --- a/examples/reactmicp/unsaturated/acc_carbo.cpp +++ b/examples/reactmicp/unsaturated/acc_carbo.cpp @@ -1,588 +1,588 @@ /* ============================================================================= 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 acc_carbo.cpp //! \brief Accelerated carbonation example //! //! This is an example to show the features of the unsaturated //! reactive transport solver //! It is a simple study of the accelerated carbonation of cement paste //! (with drying) //! //! #include "reactmicp/reactmicp_unsaturated.hpp" #include "reactmicp/systems/unsaturated/variables.hpp" #include "reactmicp/io/hdf5_unsaturated.hpp" #include "dfpm/io/hdf5_mesh.hpp" #include "specmicp_database/io/hdf5_database.hpp" #include "dfpm/solver/parabolic_structs.hpp" #include "specmicp_common/physics/unsaturated_laws.hpp" #include "specmicp_common/cli/parser.hpp" #include #include using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems::unsaturated; using VariablesBase = solver::VariablesBase; using StaggerReturnCode = solver::StaggerReturnCode; class UpscalingStagger; database::RawDatabasePtr get_database(); AdimensionalSystemSolution get_initial_condition( database::RawDatabasePtr the_database, const units::UnitsSet& units_set ); void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species); // Upscaling stagger -> user models class UpscalingStagger: public solver::UpscalingStaggerBase { public: UpscalingStagger() {} ~UpscalingStagger() {} //! \brief Initialize the stagger at the beginning of the computation //! //! \param var a shared_ptr to the variables virtual void initialize(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, VariablesBase * const var) override; //! \brief Solve the equation for the timestep //! //! \param var a shared_ptr to the variables virtual StaggerReturnCode restart_timestep( VariablesBase * const var) override; // User models scalar_t capillary_pressure(index_t node, scalar_t saturation); scalar_t vapor_pressure(index_t node, scalar_t saturation); scalar_t relative_liquid_diffusivity(index_t node, scalar_t saturation); scalar_t relative_liquid_permeability(index_t node, scalar_t saturation); scalar_t relative_gas_diffusivity(index_t node, scalar_t saturation); scalar_t intrinsic_permeability(scalar_t porosity); scalar_t intrinsic_liquid_diffusivity(scalar_t porosity); scalar_t resistance_gas_diffusivity(scalar_t porosity); private: scalar_t m_a {37.5479e6}; scalar_t m_b {2.1684}; scalar_t m_alpha {7.27}; scalar_t m_beta {0.440}; scalar_t m_pv_sat {3170}; scalar_t m_exp_r_l_d {3.00}; scalar_t m_exp_i_g_d {2.74}; scalar_t m_exp_r_g_d {4.20}; scalar_t m_exp_r_l_k {0.4191}; scalar_t m_k_0 {1e-21*1e4}; scalar_t m_d_0 {2.3e-13*1e4}; scalar_t m_phi_0 {0.2}; }; // implementations int main(int argc, char* argv[]) { cli::CommandLineParser cli_parser; cli_parser.add_option('n', "min_dt", 0.01, "minimum timestep (s)"); cli_parser.add_option('m', "max_dt", 10.0, "maximum timestep (s)"); cli_parser.add_option('d', "dx", 0.05, "space step (cm)"); cli_parser.add_option('r', "run_until", 10.0*24.0*3600.0, "run simulation until s"); cli_parser.add_option('p', "p_h20", 1000.0, "vapor pressure at the boundary"); cli_parser.add_option('c', "p_co2", 500.0, "vapor pressure at the boundary"); cli_parser.set_help_message("Atmospheric carbonation"); cli_parser.parse(argc, argv); init_logger(&std::cout, logger::Warning); index_t nb_nodes = 25; scalar_t dx = cli_parser.get_option("dx"); // cm scalar_t cross_section = 5; // cm^2 units::UnitsSet units_set; units_set.length = units::LengthUnit::centimeter; specmicp::RawDatabasePtr the_database = get_database(); AdimensionalSystemSolution init = get_initial_condition( the_database, units_set); AdimensionalSystemSolutionExtractor extr(init, the_database, units_set); std::cout << extr.porosity() << " - " << extr.saturation_water() << std::endl; //return EXIT_SUCCESS; mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d({dx, nb_nodes, cross_section}); index_t id_co2 = the_database->get_id_component("CO2"); VariablesInterface variable_interface(the_mesh, the_database, {0, id_co2}, units_set); variable_interface.set_binary_gas_diffusivity(0, 0.240); variable_interface.set_binary_gas_diffusivity(id_co2, 0.160); std::shared_ptr upscaling_stagger = std::make_shared(); auto* call_ptr = upscaling_stagger.get(); variable_interface.set_vapor_pressure_model( [call_ptr](index_t node, scalar_t sat) -> scalar_t { return call_ptr->vapor_pressure(node, sat); }); for (index_t node=1; nodenb_nodes(); ++node) { variable_interface.initialize_variables(node, extr); } variable_interface.set_porosity(0, 1.0); variable_interface.set_liquid_saturation(0, 0.0); variable_interface.set_partial_pressure(id_co2, 0, cli_parser.get_option("p_co2")); variable_interface.set_partial_pressure(0, 0, cli_parser.get_option("p_h20")); std::shared_ptr vars = variable_interface.get_variables(); vars->set_aqueous_scaling(0, 1.0e3/18.3e-3*1e-6); for (auto aqc: the_database->range_aqueous_component()) { vars->set_aqueous_scaling(aqc, 100e-6); } //vars->set_gaseous_scaling(0, 1.0/vars->get_rt()); //vars->set_gaseous_scaling(id_co2, 1.0/vars->get_rt()); std::cout << "Youpla : " << vars->get_rt() << std::endl; auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), the_database->nb_component()); bcs->add_gas_node(0); bcs->fork_constraint(0, "gas_node"); auto chem_options = EquilibriumOptionsVector::make(the_mesh->nb_nodes()); auto& def_chem_opts = chem_options->get("default"); def_chem_opts.units_set = units_set; micpsolver::MiCPSolverOptions& copts = def_chem_opts.solver_options; copts.set_maximum_iterations(200); copts.set_maximum_step_length(10, 200); copts.set_tolerance(1e-8, 1e-12); def_chem_opts.system_options.non_ideality_tolerance = 1e-14; std::shared_ptr chemistry_stagger = EquilibriumStagger::make(vars, bcs, chem_options); upscaling_stagger->initialize(vars.get()); vars->set_relative_variables(0); chemistry_stagger->initialize(vars.get()); std::shared_ptr transport_stagger = UnsaturatedTransportStagger::make(vars, bcs, true, true); transport_stagger->initialize(vars.get()); bcs->get_constraint("gas_node").set_water_partial_pressure_model( std::bind(&UpscalingStagger::vapor_pressure, upscaling_stagger.get(), 1, std::placeholders::_1) ); auto check = variable_interface.check_variables(); if (check >= VariablesValidity::error) { std::cout << variable_interface.get_liquid_saturation() << std::endl; throw std::runtime_error("Invalid variables"); } solver::ReactiveTransportSolver react_solver(transport_stagger, chemistry_stagger, upscaling_stagger); solver::ReactiveTransportOptions& opts = react_solver.get_options(); opts.residuals_tolerance = 1e-2; opts.good_enough_tolerance = 0.99; opts.maximum_iterations = 100; transport_stagger->get_saturation_options()->maximum_iterations = 500; transport_stagger->get_saturation_options()->step_tolerance = 1e-10; transport_stagger->get_saturation_options()->residuals_tolerance = 1e-8; for (auto ind: the_database->range_aqueous_component()) { auto* opts = transport_stagger->get_aqueous_options(ind); opts->step_tolerance = 1e-10; opts->residuals_tolerance = 1e-8; } //transport_stagger->get_gas_options(3)->maximum_iterations = 5000; opts.step_tolerance = 1e-10; io::UnsaturatedHDF5Saver saver("/tmp/out_acc_carbo.hdf5", vars, units_set); saver.save_timestep(0.0); solver::SimulationInformation simul_info("acc_carbo", 500); auto out_pol = [&saver]( scalar_t timestep, solver::VariablesBasePtr _) { saver.save_timestep(timestep); }; database::Database db_handler(the_database); db_handler.save("acc_carbo_data.yml"); //solver::ReactiveTransportReturnCode retcode = react_solver.solve_timestep(0.01, vars); //retcode = react_solver.solve_timestep(1.0, vars); //std::cout << (int) retcode << std::endl; // retcode = react_solver.solve_timestep(5.0, vars); // std::cout << (int) retcode << std::endl; // retcode = react_solver.solve_timestep(10.0, vars); // std::cout << (int) retcode << std::endl; // retcode = react_solver.solve_timestep(10.0, vars); // std::cout << (int) retcode << std::endl; solver::ReactiveTransportRunner runner( react_solver, cli_parser.get_option("min_dt"), cli_parser.get_option("max_dt"), simul_info); runner.set_output_policy(out_pol); runner.run_until(cli_parser.get_option("run_until"), vars); std::cout << vars->get_liquid_saturation().variable << std::endl; std::cout << vars->get_liquid_saturation()(1) - vars->get_liquid_saturation()(2) << std::endl; std::cout << vars->get_pressure_main_variables(id_co2).variable << std::endl; } void UpscalingStagger::initialize(VariablesBase * const var) { UnsaturatedVariables * const true_var = static_cast(var); // set relative models true_var->set_capillary_pressure_model( [this](index_t x, scalar_t sat){ return this->capillary_pressure(x, sat); }); true_var->set_relative_liquid_permeability_model( [this](index_t x, scalar_t sat){ return this->relative_liquid_permeability(x, sat); }); true_var->set_relative_liquid_diffusivity_model( [this](index_t x, scalar_t sat){ return this->relative_liquid_diffusivity(x, sat); }); true_var->set_relative_gas_diffusivity_model( [this](index_t x, scalar_t sat){ return this->relative_gas_diffusivity(x, sat); }); true_var->set_vapor_pressure_model( [this](index_t x, scalar_t sat){ return this->vapor_pressure(x, sat); }); // initialization SecondaryTransientVariable& porosity = true_var->get_porosity(); SecondaryVariable& permeability = true_var->get_liquid_permeability(); SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity(); SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity(); true_var->set_relative_variables(); gas_diffusivity(0) = resistance_gas_diffusivity(1.0); true_var->get_relative_gas_diffusivity()(0) = relative_gas_diffusivity(0, 1.0); for (index_t node=1; nodeget_mesh()->nb_nodes(); ++node) { scalar_t phi = porosity(node); permeability(node) = intrinsic_permeability(phi); liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi); gas_diffusivity(node) = resistance_gas_diffusivity(phi); //true_var->set_relative_variables(node); } } scalar_t UpscalingStagger::capillary_pressure(index_t node, scalar_t saturation) { if (saturation == 0.0) return 0.0; return laws::capillary_pressure_BaroghelBouny(saturation, m_a, m_b); } scalar_t UpscalingStagger::vapor_pressure(index_t node, scalar_t saturation) { scalar_t pv; if (saturation == 0.0) pv = 0.0; else if (saturation >= 1.0) pv = m_pv_sat; else pv = m_pv_sat*laws::invert_kelvin_equation( capillary_pressure(node, saturation)); return pv; } scalar_t UpscalingStagger::relative_liquid_diffusivity( index_t node, scalar_t saturation ) { if (saturation == 0.0) return 0.0; return std::pow(saturation, m_exp_r_l_d); } scalar_t UpscalingStagger::relative_liquid_permeability( index_t node, scalar_t saturation ) { if (saturation == 0.0) return 0.0; return laws::relative_liquid_permeability_Mualem(saturation, 1.0/m_b); } scalar_t UpscalingStagger::relative_gas_diffusivity(index_t node, scalar_t saturation) { if (saturation == 0.0) return 1.0; return std::pow(1.0 - saturation, m_exp_r_g_d); } scalar_t UpscalingStagger::intrinsic_permeability(scalar_t porosity) { scalar_t tmp_1 = porosity / m_phi_0; scalar_t tmp_2 = (1.0 - m_phi_0) / (1.0 - porosity); tmp_1 = std::pow(tmp_1, 3); tmp_2 = std::pow(tmp_2, 2); return m_k_0 * tmp_1 * tmp_2; } scalar_t UpscalingStagger::intrinsic_liquid_diffusivity(scalar_t porosity) { return m_d_0 * std::exp(-9.95*porosity); } scalar_t UpscalingStagger::resistance_gas_diffusivity(scalar_t porosity) { return std::pow(porosity, m_exp_i_g_d); } void UpscalingStagger::initialize_timestep( scalar_t dt, VariablesBase * const var ) { } StaggerReturnCode UpscalingStagger::restart_timestep(VariablesBase * const var) { UnsaturatedVariables* true_var = static_cast(var); SecondaryTransientVariable& porosity = true_var->get_porosity(); SecondaryVariable& permeability = true_var->get_liquid_permeability(); SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity(); SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity(); for (index_t node=1; nodeget_mesh()->nb_nodes(); ++node) { scalar_t phi = porosity(node); permeability(node) = intrinsic_permeability(phi); liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi); gas_diffusivity(node) = resistance_gas_diffusivity(phi); true_var->set_relative_variables(node); } return StaggerReturnCode::ResidualMinimized; } specmicp::RawDatabasePtr get_database() { specmicp::database::Database thedatabase(CEMDATA_PATH); std::map swap = {{"HCO3[-]", "CO2"},}; thedatabase.swap_components(swap); thedatabase.remove_half_cell_reactions(); database::RawDatabasePtr raw_data = thedatabase.get_database(); return raw_data; } AdimensionalSystemSolution get_initial_condition( database::RawDatabasePtr the_database, const units::UnitsSet& units_set ) { Vector oxyde_compo(4); oxyde_compo << 67.98, 22.66, 5.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 = 5e-4; 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*0.0298}, {"KOH", mult*0.0801}, {"Cl[-]", mult*0.0001} }; Dissolver dissolver(the_database); UpscalingStagger upsstag_init; // just for vapor pressure model AdimensionalSystemConstraints constraints; constraints.set_total_concentrations( dissolver.dissolve(formulation, true) ); //constraints.set_inert_volume_fraction(0.1); AdimensionalSystemSolver the_solver(the_database, constraints); the_solver.get_options().units_set = units_set; the_solver.get_options().system_options.non_ideality_tolerance = 1e-12; micpsolver::MiCPSolverOptions* solver_options = &the_solver.get_options().solver_options; //solver_options.maxstep = 10.0; solver_options->maxiter_maxstep = 200; solver_options->max_iter = 200; solver_options->set_tolerance(1e-9); Vector x; - the_solver.initialise_variables(x, 0.8, -4); + the_solver.initialize_variables(x, 0.8, -4); micpsolver::MiCPPerformance perf = the_solver.solve(x); if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success) { throw std::runtime_error("Error : failed to solve initial conditions"); } auto prev_sol = the_solver.get_raw_solution(x); prev_sol.main_variables(0) *= 0.8; constraints.set_water_partial_pressure_model( std::bind(&UpscalingStagger::vapor_pressure, &upsstag_init, 1, std::placeholders::_1) ); the_solver = AdimensionalSystemSolver (the_database, constraints, prev_sol); the_solver.get_options().units_set = units_set; the_solver.get_options().system_options.non_ideality_tolerance = 1e-12; solver_options = &the_solver.get_options().solver_options; //solver_options.maxstep = 10.0; solver_options->maxiter_maxstep = 200; solver_options->set_tolerance(1e-9); perf = the_solver.solve(x); if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success) { throw std::runtime_error("Error : failed to solve initial conditions"); } return the_solver.get_raw_solution(x); } 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); } diff --git a/examples/specmicp/adimensional/carbofe.cpp b/examples/specmicp/adimensional/carbofe.cpp index 82f3e28..9af9447 100644 --- a/examples/specmicp/adimensional/carbofe.cpp +++ b/examples/specmicp/adimensional/carbofe.cpp @@ -1,223 +1,223 @@ /* ============================================================================= 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 "specmicp/specmicp.hpp" #include "specmicp_common/timer.hpp" #include class AutomaticReactionPath: public specmicp::EquilibriumCurve { public: AutomaticReactionPath() { specmicp::database::Database thedatabase("../data/cemdata.yaml"); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"}, {"Al[3+]","Al(OH)4[-]"}, {"Fe[2+]", "Fe(OH)4[-]"}, }); thedatabase.swap_components(swapping); thedatabase.remove_half_cell_reactions(std::vector({"SO4[2-]"})); thedatabase.remove_gas_phases(); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); set_database(raw_data); specmicp::Formulation formulation; specmicp::scalar_t mult = 6e3; specmicp::scalar_t m_c3s = mult*0.6; specmicp::scalar_t m_c2s = mult*0.2; specmicp::scalar_t m_c3a = mult*0.10; specmicp::scalar_t m_c4af = mult*0.05; specmicp::scalar_t m_gypsum = mult*0.05; specmicp::scalar_t wc = 0.8; specmicp::scalar_t m_water = wc*1e-3*( m_c3s*(3*56.08+60.08) + m_c2s*(2*56.06+60.08) + m_c3a*(3*56.08+101.96) + m_c4af*(4*56.08+101.96+159.07) + m_gypsum*(56.08+80.06+2*18.02) ); formulation.mass_solution = m_water; formulation.amount_minerals = { {"C3S", m_c3s}, {"C2S", m_c2s}, {"C3A", m_c3a}, {"C4AF", m_c4af}, {"Gypsum", m_gypsum} }; formulation.extra_components_to_keep = {"HCO3[-]", }; formulation.minerals_to_keep = { "Portlandite", "CSH,jennite", "CSH,tobermorite", "SiO2,am", "Calcite", "Al(OH)3,am", "C3AH6", "C4AH13", "C2AH8", "CAH10", "Monosulfoaluminate", "Tricarboaluminate", "Monocarboaluminate", "Hemicarboaluminate", //"Straetlingite", "Gypsum", "Ettringite", "Thaumasite", "C3FH6", "C4FH13", "C2FH8", "Fe(OH)3(mic)", "Fe-Ettringite", "Fe-Monosulfate", "Fe-Monocarbonate", "Fe-Hemicarbonate" }; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); id_h2o = thedatabase.component_label_to_id("H2O"); id_ho = thedatabase.component_label_to_id("HO[-]"); id_hco3 = thedatabase.component_label_to_id("HCO3[-]"); id_co2g = thedatabase.gas_label_to_id("CO2(g)"); id_ca = thedatabase.component_label_to_id("Ca[2+]"); constraints() = specmicp::AdimensionalSystemConstraints(total_concentrations); constraints().charge_keeper = id_ho; specmicp::AdimensionalSystemSolverOptions& options = solver_options(); options.solver_options.maxstep = 50.0; options.solver_options.max_iter = 200; options.solver_options.maxiter_maxstep = 200; options.solver_options.use_crashing = false; options.solver_options.use_scaling = true; options.solver_options.non_monotone_linesearch = true; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 300; options.solver_options.projection_min_variable = 1e-9; options.solver_options.fvectol = 1e-10; options.solver_options.steptol = 1e-16; options.system_options.non_ideality_tolerance = 1e-14; options.system_options.scaling_electron = 1e10; specmicp::Vector x; specmicp::AdimensionalSystemSolver solver(raw_data, constraints(), solver_options()); - solver.initialise_variables(x, 0.5, { + solver.initialize_variables(x, 0.5, { {"Ca[2+]", -2.0}, {"HO[-]", -1.4}, {"Al(OH)4[-]", -4}, {"Fe(OH)4[-]", -4}, {"SiO(OH)3[-]", -6.0}, {"SO4[2-]", -6.2} }); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); if (perf.return_code <= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet) { error_handling("Failed to solve first problem, return code : " + std::to_string((int) perf.return_code) + "."); } solution_vector() = x; initialize_solution(solver.get_raw_solution(x)); // std::cout << x << std::endl; // std::cout << current_solution().secondary_molalities << std::endl; std::cout << "Buffering \t" << "porosity \t pH \t Eh"; for (specmicp::index_t mineral: raw_data->range_mineral()) { std::cout << "\t" << raw_data->minerals.get_label(mineral); } std::cout << std::endl; output(); } void output() { specmicp::AdimensionalSystemSolutionExtractor sol(current_solution(), database(), solver_options().units_set); std::cout << constraints().total_concentrations(id_hco3)/constraints().total_concentrations(id_ca) << "\t" << sol.porosity() << "\t" << sol.pH() << "\t" << sol.Eh(); for (specmicp::index_t mineral: database()->range_mineral()) { std::cout << "\t" << sol.volume_fraction_mineral(mineral); } std::cout << std::endl; } void update_problem() { constraints().total_concentrations(id_hco3) += 100; constraints().total_concentrations(id_ho) -= 100; constraints().total_concentrations(id_h2o) += 100; } private: specmicp::index_t id_h2o; specmicp::index_t id_ho; specmicp::index_t id_hco3; specmicp::index_t id_co2g; specmicp::index_t id_ca; }; int main() { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; specmicp::Timer timer; AutomaticReactionPath test_automatic; for (int i=0; i<180; ++i) { test_automatic.run_step(); } timer.stop(); std::cout << "Execution time : " << timer.elapsed_time() << "s" << std::endl; return EXIT_SUCCESS; } diff --git a/examples/specmicp/adimensional/momas_thermo.cpp b/examples/specmicp/adimensional/momas_thermo.cpp index f271367..0194ead 100644 --- a/examples/specmicp/adimensional/momas_thermo.cpp +++ b/examples/specmicp/adimensional/momas_thermo.cpp @@ -1,284 +1,284 @@ /* ============================================================================= 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 "specmicp_common/log.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp/adimensional/equilibrium_curve.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional_kinetics/kinetic_variables.hpp" #include "specmicp/adimensional_kinetics/kinetic_model.hpp" #include "specmicp/adimensional_kinetics/kinetic_system_euler_solver.hpp" #include "specmicp/adimensional_kinetics/kinetic_system_solver.hpp" #include "specmicp_database/database.hpp" #include "specmicp_database/aqueous_selector.hpp" #include specmicp::RawDatabasePtr get_momas_db() { specmicp::database::Database thedatabase("../data/momas_benchmark.yaml"); thedatabase.remove_components({"X5",}); specmicp::database::AqueousSelector selector(thedatabase.get_database()); selector.remove_aqueous({selector.aqueous_label_to_id("C6"), selector.aqueous_label_to_id("C7")}); thedatabase.save("out_momas.db"); return thedatabase.get_database(); } class AutomaticReactionPath: public specmicp::EquilibriumCurve { public: AutomaticReactionPath() { specmicp::RawDatabasePtr raw_data = get_momas_db(); set_database(raw_data); specmicp::Vector total_concentrations(raw_data->nb_component()); total_concentrations << 1, 0.0, -3.0, 0.0, 1.0; constraints() = specmicp::AdimensionalSystemConstraints(total_concentrations); constraints().disable_conservation_water(); constraints().surface_model.model_type = specmicp::SurfaceEquationType::Equilibrium; constraints().surface_model.concentration = 1.0; specmicp::AdimensionalSystemSolverOptions& options = solver_options(); options.solver_options.maxstep = 1.0; options.solver_options.max_iter = 200; options.solver_options.maxiter_maxstep = 200; options.solver_options.use_crashing = false; options.solver_options.use_scaling = true; options.solver_options.non_monotone_linesearch = false; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 100; options.solver_options.projection_min_variable = 1e-9; options.solver_options.fvectol = 1e-6; options.solver_options.steptol = 1e-14; options.system_options.non_ideality_tolerance = 1e-8; options.system_options.non_ideality = false; options.use_pcfm = true; //solve_first_problem(); specmicp::Vector variables; specmicp::AdimensionalSystemSolver solver(raw_data, constraints(), solver_options()); - solver.initialise_variables(variables, + solver.initialize_variables(variables, 0.8, {{"X2", -3}, {"X4", -11}}, {}, 0.5 ); solver.run_pcfm(variables); specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, false); std::cout << (int) current_perf.return_code << std::endl; solution_vector() = variables; initialize_solution(solver.get_raw_solution(variables)); std::cout << variables << std::endl; std::cout << "Buffering \t"; for (specmicp::index_t mineral: raw_data->range_mineral()) { std::cout << "\t" << raw_data->get_label_mineral(mineral); } std::cout << "\t S"; for (specmicp::index_t sorbed: raw_data->range_sorbed()) { std::cout << "\t" << raw_data->get_label_sorbed(sorbed); } std::cout << std::endl; output(); } void output() { specmicp::AdimensionalSystemSolutionExtractor sol(current_solution(), database(), solver_options().units_set); std::cout << constraints().total_concentrations(1); for (specmicp::index_t mineral: database()->range_mineral()) { std::cout << "\t" << sol.mole_concentration_mineral(mineral); } std::cout << "\t" << sol.free_surface_concentration(); for (specmicp::index_t sorbed: database()->range_sorbed()) { std::cout << "\t" << sol.molality_sorbed_species(sorbed); } std::cout << std::endl; } void update_problem() { constraints().total_concentrations(1) += 0.1; constraints().total_concentrations(2) += 0.1; constraints().total_concentrations(3) += 0.1; } private: specmicp::index_t id_h2o; specmicp::index_t id_ho; specmicp::index_t id_hco3; specmicp::index_t id_co2g; specmicp::index_t id_ca; }; class MomasKinetics: public specmicp::kinetics::AdimKineticModel { public: MomasKinetics(): specmicp::kinetics::AdimKineticModel({0}), m_id_cc(0), m_id_c3(2), m_id_x4(4) {} //! \brief Compute the kinetic rates and store them in dydt void compute_rate( specmicp::scalar_t t, const specmicp::Vector& y, specmicp::kinetics::AdimKineticVariables& variables, specmicp::Vector& dydt ) { dydt.resizeLike(y); specmicp::AdimensionalSystemSolution& solution = variables.equilibrium_solution(); specmicp::scalar_t factor = std::pow(1e3*0.5*solution.secondary_molalities(m_id_c3),3)/ std::pow(1e3*0.5*pow10(solution.main_variables(m_id_x4)), 2); factor = 0.2*factor-1.0; specmicp::scalar_t k_constant = 10.0; if (factor >= 0) k_constant = 1e-2; dydt(0) = k_constant*factor; } private: specmicp::index_t m_id_cc; specmicp::index_t m_id_c3; specmicp::index_t m_id_x4; }; void solve_kinetics_problem() { specmicp::RawDatabasePtr raw_data = get_momas_db(); specmicp::Vector total_concentrations(6); total_concentrations << 1, 0.6, -2.4, 0.6, 1.0, 1.0; specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.disable_conservation_water(); specmicp::AdimensionalSystemSolverOptions options; options.solver_options.maxstep = 10.0; options.solver_options.max_iter = 200; options.solver_options.maxiter_maxstep = 200; options.solver_options.use_crashing = false; options.solver_options.use_scaling = true; options.solver_options.non_monotone_linesearch = true; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 50; options.solver_options.projection_min_variable = 1e-9; options.solver_options.fvectol = 1e-8; options.solver_options.steptol = 1e-10; options.system_options.non_ideality_tolerance = 1e-10; options.system_options.non_ideality = false; specmicp::Vector equil_variables; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); - solver.initialise_variables(equil_variables, 0.8, -3); + solver.initialize_variables(equil_variables, 0.8, -3); specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(equil_variables, false); std::cout << (int) current_perf.return_code << std::endl; specmicp::AdimensionalSystemSolution equil_solution = solver.get_raw_solution(equil_variables); std::shared_ptr model = std::make_shared(); specmicp::Vector mineral_concentrations(1); mineral_concentrations << 5.0; specmicp::kinetics::AdimKineticSystemEulerSolver kinetic_solver( //specmicp::kinetics::AdimKineticSystemSolver kinetic_solver( model, total_concentrations, mineral_concentrations, constraints, equil_solution, raw_data); options.solver_options.use_scaling = false; kinetic_solver.get_options().speciation_options = options; kinetic_solver.solve(0.005, 1.0); std::cout << "Final amount : " << kinetic_solver.variables().concentration_mineral(0) << std::endl; } int main() { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Debug; // solve_kinetics_problem(); AutomaticReactionPath test_automatic; for (int i=0; i<20; ++i) { test_automatic.run_step(); } return EXIT_SUCCESS; } diff --git a/src/specmicp/adimensional/adimensional_system_solver.cpp b/src/specmicp/adimensional/adimensional_system_solver.cpp index 8e96840..bb55d01 100644 --- a/src/specmicp/adimensional/adimensional_system_solver.cpp +++ b/src/specmicp/adimensional/adimensional_system_solver.cpp @@ -1,448 +1,491 @@ /* ============================================================================= 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 "adimensional_system_solver.hpp" #include "adimensional_system.hpp" #include "adimensional_system_solution.hpp" #include "adimensional_system_pcfm.hpp" #include "specmicp_common/micpsolver/micpsolver.hpp" #include "specmicp_common/log.hpp" #include namespace specmicp { AdimensionalSystemSolution solve_equilibrium( std::shared_ptr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ) { AdimensionalSystemSolver solver(data, constraints, options); Vector variables; micpsolver::MiCPPerformance perf = solver.solve(variables, true); if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet) throw std::runtime_error("Failed to solve the problem"); return solver.get_raw_solution(variables); } AdimensionalSystemSolution SPECMICP_DLL_PUBLIC solve_equilibrium( RawDatabasePtr data_ptr, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemSolverOptions& options ) { AdimensionalSystemSolver solver(data_ptr, constraints, previous_solution, options); Vector variables; micpsolver::MiCPPerformance perf = solver.solve(variables, true); if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet) throw std::runtime_error("Failed to solve the problem"); return solver.get_raw_solution(variables); } // constructor // ----------- AdimensionalSystemSolver::AdimensionalSystemSolver( RawDatabasePtr data, const AdimensionalSystemConstraints& constraints ): OptionsHandler(), m_data(data), m_system(std::make_shared( data, constraints, get_options().system_options, get_options().units_set )), m_var(Vector::Zero(data->nb_component()+1+data->nb_mineral())) {} AdimensionalSystemSolver::AdimensionalSystemSolver( RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ): OptionsHandler(options), m_data(data), m_system(std::make_shared( data, constraints, options.system_options, options.units_set )), m_var(Vector::Zero(data->nb_component()+1+data->nb_mineral())) {} AdimensionalSystemSolver::AdimensionalSystemSolver( RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution ): OptionsHandler(), m_data(data), m_system(std::make_shared( data, constraints, previous_solution, get_options().system_options, get_options().units_set )), m_var(Vector::Zero(data->nb_component()+1+data->nb_mineral())) {} AdimensionalSystemSolver::AdimensionalSystemSolver( RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemSolverOptions& options ): OptionsHandler(options), m_data(data), m_system(std::make_shared( data, constraints, previous_solution, options.system_options, options.units_set )), m_var(Vector::Zero(data->nb_component()+1+data->nb_mineral())) {} AdimensionalSystemSolution AdimensionalSystemSolver::get_raw_solution(Vector& x) { set_true_variable_vector(x); return m_system->get_solution(x, m_var); } // Solving the system // ------------------ micpsolver::MiCPPerformance AdimensionalSystemSolver::solve(Vector& x, bool init) { m_system->set_units(get_options().units_set); if (init) { m_system->reasonable_starting_guess(x); if (get_options().use_pcfm) { run_pcfm(x); } } else if (get_options().force_pcfm) { run_pcfm(x); } set_true_variable_vector(x); micpsolver::MiCPPerformance perf = solve_system(); int cnt = 0; while (perf.return_code < micpsolver::MiCPSolverReturnCode::Success // != micpsolver::MiCPSolverReturnCode::ResidualMinimized and get_options().allow_restart and cnt < 3) { WARNING << "Failed to solve the system ! Return code :" << (int) perf.return_code << ". We shake it up and start again"; const scalar_t save_penalization_factor = get_options().solver_options.penalization_factor; const bool save_scaling = get_options().solver_options.use_scaling; get_options().solver_options.use_scaling = true; if (save_penalization_factor == 1) get_options().solver_options.penalization_factor = 0.8; set_return_vector(x); m_system->reasonable_restarting_guess(x); if (get_options().use_pcfm or get_options().force_pcfm) run_pcfm(x); set_true_variable_vector(x); micpsolver::MiCPPerformance perf2 = solve_system(); get_options().solver_options.penalization_factor = save_penalization_factor; get_options().solver_options.use_scaling = save_scaling; perf += perf2; ++cnt; } if (perf.return_code > micpsolver::MiCPSolverReturnCode::NotConvergedYet) set_return_vector(x); return perf; } micpsolver::MiCPPerformance AdimensionalSystemSolver::solve_system() { micpsolver::MiCPSolver solver(m_system); solver.set_options(get_options().solver_options); solver.solve(m_var); return solver.get_perfs(); } // Variables management // --------------------- void AdimensionalSystemSolver::set_true_variable_vector(const Vector& x) { const std::vector& non_active_component = m_system->get_non_active_component(); uindex_t new_i {0}; if (m_system->is_active_component(0)) { m_var(new_i) = x(m_system->dof_water()); ++new_i; } for (index_t i: m_data->range_aqueous_component()) { const auto it = std::find(non_active_component.cbegin(), non_active_component.cend(),i); if (it != non_active_component.cend()) continue; scalar_t value = x(m_system->dof_component(i)); if (value == -HUGE_VAL) // check for previously undefined value { value = m_system->get_options().new_component_concentration; } m_var(new_i) = value; ++new_i; } if (m_system->ideq_surf() != no_equation) { m_var(new_i) = x(m_system->dof_surface()); ++new_i; } if (m_system->ideq_electron() != no_equation) { m_var(new_i) = x(m_system->dof_electron()); ++new_i; } for (index_t m: m_data->range_mineral()) { bool to_keep = true; for (auto it=non_active_component.begin(); it!=non_active_component.end(); ++it) { if (m_data->nu_mineral(m, *it) != 0) to_keep = false; } if (to_keep) { m_var(new_i) = x(m_system->dof_mineral(m)); ++new_i; } } m_var.conservativeResize(new_i); specmicp_assert(new_i == (unsigned) m_system->total_variables()); } void AdimensionalSystemSolver::set_return_vector(Vector& x) { const std::vector& non_active_component = m_system->get_non_active_component(); uindex_t new_i = 0; if (m_system->is_active_component(0)) { x(m_system->dof_water()) = m_var(new_i); ++new_i; } else { x(m_system->dof_water()) = m_system->volume_fraction_water(x); } for (index_t i: m_data->range_aqueous_component()) { const auto it = std::find(non_active_component.cbegin(), non_active_component.cend(),i); if (it != non_active_component.cend()) { x(m_system->dof_component(i)) = -HUGE_VAL; continue; } x(m_system->dof_component(i)) = m_var(new_i) ; ++new_i; } if (m_system->ideq_surf() != no_equation) { x(m_system->dof_surface()) = m_var(new_i); ++new_i; } else x(m_system->dof_surface()) = -HUGE_VAL; if (m_system->ideq_electron() != no_equation) { x(m_system->dof_electron()) = m_var(new_i); ++new_i; } else x(m_system->dof_electron()) = -HUGE_VAL; for (index_t m: m_data->range_mineral()) { bool to_keep = true; for (const index_t& k: non_active_component) { if (m_data->nu_mineral(m, k) != 0.0) to_keep = false; } if (to_keep) { x(m_system->dof_mineral(m)) =m_var(new_i); ++new_i; } else { x(m_system->dof_mineral(m)) = 0.0; } } } // PCFM // ---- void AdimensionalSystemSolver::run_pcfm(Vector &x) { DEBUG << "Start PCFM initialization."; // we set up the true variable set_true_variable_vector(x); // The residual is computed to have some point of comparison Vector residuals(m_system->total_variables()); residuals.setZero(); m_system->get_residuals(m_var, residuals); const scalar_t res_0 = residuals.norm(); // the pcfm iterations are executed AdimensionalSystemPCFM pcfm_solver(m_system); PCFMReturnCode retcode = pcfm_solver.solve(m_var); // Check the answer if (retcode < PCFMReturnCode::Success) { // small prograss is still good enough m_system->get_residuals(m_var, residuals); const scalar_t final_residual = residuals.norm(); DEBUG << "Final pcfm residuals : " << final_residual << " set_secondary_variables(m_var); } } // Initialisation of variables // --------------------------- void AdimensionalSystemSolver::initialise_variables( Vector& x, scalar_t volume_fraction_water, std::map log_molalities, std::map volume_fraction_minerals, scalar_t log_free_sorption_site_concentration ) { + std::unordered_map un_mol( + log_molalities.begin(), log_molalities.end()); + std::unordered_map un_sol( + volume_fraction_minerals.begin(), + volume_fraction_minerals.end() + ); + initialize_variables(x, + volume_fraction_water, + un_mol, un_sol, + log_free_sorption_site_concentration + ); +} + +void AdimensionalSystemSolver::initialize_variables( + Vector& x, + scalar_t volume_fraction_water, + std::unordered_map log_molalities, + std::unordered_map volume_fraction_minerals, + scalar_t log_free_sorption_site_concentration + ) +{ + m_system->reasonable_starting_guess(x); if (volume_fraction_water < 0 or volume_fraction_water > 1) { WARNING << "Initial guess for the volume fraction of water is not between 0 and 1"; } x(m_system->dof_water()) = volume_fraction_water; for (auto pair: log_molalities) { index_t idc = m_data->get_id_component(pair.first); if (idc == no_species or idc == m_data->electron_index() or idc == m_data->water_index()) { throw std::invalid_argument("This is not an aqueous component : "+pair.first); } if (pair.second > 0) { WARNING << "Initial molality for : " << pair.first << "is bigger than 1 molal."; } x(m_system->dof_component(idc)) = pair.second; } for (auto pair: volume_fraction_minerals) { index_t idm = m_data->get_id_mineral(pair.first); if (idm == no_species ) { throw std::invalid_argument("This is not a mineral at equilibrium : "+pair.first); } if (pair.second < 0 or pair.second > 1) { WARNING << "Initial volume fraction for : " << pair.first << "is not between 0 and 1."; } x(m_system->dof_mineral(idm)) = pair.second; } if (log_free_sorption_site_concentration != 0.0) x(m_system->dof_surface()) = log_free_sorption_site_concentration; +} +void AdimensionalSystemSolver::initialise_variables( + Vector& x, + scalar_t volume_fraction_water, + scalar_t log_molalities + ) +{ + initialize_variables(x, volume_fraction_water, log_molalities); } -void AdimensionalSystemSolver::initialise_variables(Vector& x, - scalar_t volume_fraction_water, - scalar_t log_molalities - ) +void AdimensionalSystemSolver::initialize_variables( + Vector& x, + scalar_t volume_fraction_water, + scalar_t log_molalities + ) { m_system->reasonable_starting_guess(x); if (volume_fraction_water < 0 or volume_fraction_water > 1) { WARNING << "Initial guess for the volume fraction of water is not between 0 and 1"; } x(m_system->dof_water()) = volume_fraction_water; if (log_molalities > 0) { WARNING << "Initial molality for : " << log_molalities << "is bigger than 1 molal."; } x.segment(1, m_data->nb_component()-1).setConstant(log_molalities); } +void AdimensionalSystemSolver::initialise_variables( + Vector& x, + scalar_t volume_fraction_water, + scalar_t log_molalities, + scalar_t log_free_sorption_site_concentration + ) +{ + initialize_variables(x, volume_fraction_water, + log_molalities, log_free_sorption_site_concentration); +} -void AdimensionalSystemSolver::initialise_variables(Vector& x, - scalar_t volume_fraction_water, - scalar_t log_molalities, - scalar_t log_free_sorption_site_concentration - ) +void AdimensionalSystemSolver::initialize_variables( + Vector& x, + scalar_t volume_fraction_water, + scalar_t log_molalities, + scalar_t log_free_sorption_site_concentration + ) { initialise_variables(x, volume_fraction_water, log_molalities); x(m_system->dof_surface()) = log_free_sorption_site_concentration; } + } // end namespace specmicp diff --git a/src/specmicp/adimensional/adimensional_system_solver.hpp b/src/specmicp/adimensional/adimensional_system_solver.hpp index 80c051c..265e7d1 100644 --- a/src/specmicp/adimensional/adimensional_system_solver.hpp +++ b/src/specmicp/adimensional/adimensional_system_solver.hpp @@ -1,205 +1,245 @@ /* ============================================================================= 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_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP #define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP #include "adimensional_system_solver_structs.hpp" #include "specmicp_database/database.hpp" #include "specmicp_common/options_handler.hpp" +#include #include //! \file adimensional_system_solver.hpp //! \brief Solve a reduced system namespace specmicp { // forward declaration class AdimensionalSystem; struct AdimensionalSystemSolution; //! \brief Solve an adimensional system //! //! Take care of non-existing component in the system //! and restart the computation if necessary //! //! \ingroup specmicp_api class SPECMICP_DLL_PUBLIC AdimensionalSystemSolver: public OptionsHandler { public: using SystemPtr = std::shared_ptr; //! Default constructor //! //! Another constructor is necessary AdimensionalSystemSolver() {} //! \brief Initialise a solver from scratch //! //! \param data A raw database //! \param constraints The system to solver AdimensionalSystemSolver(RawDatabasePtr data, const AdimensionalSystemConstraints& constraints ); //! \brief Initialise a solver from scratch //! //! \param data A raw database //! \param constraints The system to solver //! \param options (optional) customize the behavior of the solver AdimensionalSystemSolver(RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ); //! \brief Initialise a solver from a previous solution //! //! \param data A raw database //! \param constraints The system to solver //! \param previous_solution A solution of a similar system that will be used to initialize the system AdimensionalSystemSolver(RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution ); //! \brief Initialise a solver from a previous solution //! //! \param data A raw database //! \param constraints The system to solver //! \param previous_solution A solution of a similar system that will be used to initialize the system //! \param options customize the behavior of the solver AdimensionalSystemSolver(RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemSolverOptions& options ); //! \brief solve the problem using initial guess x //! //! \param[in,out] x in -> initial guess, out -> solution //! \param init if true, the algorithm guess a starting point micpsolver::MiCPPerformance solve(Vector& x, bool init=false); //! \brief Return the system used for the computation SystemPtr get_system() {return m_system;} //! \brief Return the solution in a manageable form //! //! \param x The solution (complete set of the main variables) AdimensionalSystemSolution get_raw_solution(Vector& x); //! \brief Initialize the problem using the Positive continuous fraction method //! //! \sa specmicp::AdimensionalSystemPCFM void run_pcfm(Vector& x); //! \brief Custom initialisation of variables //! //! The amount compon //! \param[out] x The initial guess (complete set of the main variables) //! \param volume_fraction_water volume fraction of water //! \param log_molalities log_10 of the molalities for chosen aqueous component //! \param volume_fraction_minerals volume fraction of the minerals //! \param log_free_sorption_site_concentration concentration of free sorption site void initialise_variables(Vector& x, scalar_t volume_fraction_water, std::map log_molalities, std::map volume_fraction_minerals = {}, scalar_t log_free_sorption_site_concentration = 0 + ) + SPECMICP_DEPRECATED("Use initialize_variables"); + //! \brief Custom initialisation of variables + //! + //! The amount compon + //! \param[out] x The initial guess (complete set of the main variables) + //! \param volume_fraction_water volume fraction of water + //! \param log_molalities log_10 of the molalities for chosen aqueous component + //! \param volume_fraction_minerals volume fraction of the minerals + //! \param log_free_sorption_site_concentration concentration of free sorption site + void initialize_variables(Vector& x, + scalar_t volume_fraction_water, + std::unordered_map log_molalities, + std::unordered_map volume_fraction_minerals = {}, + scalar_t log_free_sorption_site_concentration = 0 ); //! \brief Custom initialisation of variables //! //! The amount compon //! \param[out] x The initial guess (complete set of the main variables) //! \param volume_fraction_water volume fraction of water //! \param log_molalities log_10 of the molalities for all aqueous components void initialise_variables(Vector& x, + scalar_t volume_fraction_water, + scalar_t log_molalities + ) + SPECMICP_DEPRECATED("Use initialize_variables"); + //! \brief Custom initialisation of variables + //! + //! The amount compon + //! \param[out] x The initial guess (complete set of the main variables) + //! \param volume_fraction_water volume fraction of water + //! \param log_molalities log_10 of the molalities for all aqueous components + void initialize_variables(Vector& x, scalar_t volume_fraction_water, scalar_t log_molalities ); //! \brief Custom initialisation of variables //! //! The amount compon //! \param[out] x The initial guess (complete set of the main variables) //! \param volume_fraction_water volume fraction of water //! \param log_molalities log_10 of the molalities for all aqueous components //! \param log_free_sorption_site_concentration concentration of free sorption site void initialise_variables(Vector& x, + scalar_t volume_fraction_water, + scalar_t log_molalities, + scalar_t log_free_sorption_site_concentration + ) + SPECMICP_DEPRECATED("Use initialize_variables"); + //! \brief Custom initialisation of variables + //! + //! The amount compon + //! \param[out] x The initial guess (complete set of the main variables) + //! \param volume_fraction_water volume fraction of water + //! \param log_molalities log_10 of the molalities for all aqueous components + //! \param log_free_sorption_site_concentration concentration of free sorption site + void initialize_variables(Vector& x, scalar_t volume_fraction_water, scalar_t log_molalities, scalar_t log_free_sorption_site_concentration ); private: //! \brief set up the true variable vector //! //! \param x The solution (complete set of the main variables) void SPECMICP_DLL_LOCAL set_true_variable_vector(const Vector& x); //! \brief set up the true solution vector //! //! add zero components //! \param x The solution (complete set of the main variables) void SPECMICP_DLL_LOCAL set_return_vector(Vector& x); //! \brief solve the problem micpsolver::MiCPPerformance SPECMICP_DLL_LOCAL solve_system(); RawDatabasePtr m_data; //! The raw database SystemPtr m_system; //! The system to solve Vector m_var; //! Copy of the solution vector (necessary in case of failing) }; //! \brief Solve a reduced system, function provided for convenience //! //! \param data_ptr the database //! \param constraints Constraints applied to the system //! \param options Options for the solver AdimensionalSystemSolution SPECMICP_DLL_PUBLIC solve_equilibrium( RawDatabasePtr data_ptr, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ); //! \brief Solve a reduced system, function provided for convenience //! //! \param data_ptr the database //! \param constraints constraints applied to the system //! \param previous_solution a previous solution //! \param options options for the solver AdimensionalSystemSolution SPECMICP_DLL_PUBLIC solve_equilibrium( RawDatabasePtr data_ptr, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemSolverOptions& options ); } // end namespace specmicp #endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP diff --git a/tests/specmicp/adim/adimensional_system_carboalu.cpp b/tests/specmicp/adim/adimensional_system_carboalu.cpp index afcac5e..a469d79 100644 --- a/tests/specmicp/adim/adimensional_system_carboalu.cpp +++ b/tests/specmicp/adim/adimensional_system_carboalu.cpp @@ -1,159 +1,159 @@ #include "catch.hpp" #include "specmicp_common/log.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "specmicp_common/timer.hpp" #include "specmicp_database/database.hpp" #include TEST_CASE("Carboalu - using adimensional system ", "[Adimensional, Carbonation, carboalu]") { std::cerr.flush(); specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; SECTION("Carboalu") { std::cout << "---------------------\n Carboalu\n ---------------------" << std::endl; specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"}, {"Al[3+]","Al(OH)4[-]"}, }); thedatabase.swap_components(swapping); //thedatabase.remove_half_cell_reactions(); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); thedatabase.remove_half_cell_reactions({"H2O", "HO[-]", "HCO3[-]"}); specmicp::Formulation formulation; specmicp::scalar_t mult = 5e3; specmicp::scalar_t m_c3s = mult*0.6; specmicp::scalar_t m_c2s = mult*0.2; specmicp::scalar_t m_c3a = mult*0.10; specmicp::scalar_t m_gypsum = mult*0.10; specmicp::scalar_t wc = 0.8; specmicp::scalar_t m_water = wc*1e-3*( m_c3s*(3*56.08+60.08) + m_c2s*(2*56.06+60.08) + m_c3a*(3*56.08+101.96) + m_gypsum*(56.08+80.06+2*18.02) ); formulation.mass_solution = m_water; formulation.amount_minerals = { {"C3S", m_c3s}, {"C2S", m_c2s}, {"C3A", m_c3a}, {"Gypsum", m_gypsum} }; formulation.extra_components_to_keep = {"HCO3[-]", }; formulation.minerals_to_keep = { "Portlandite", "CSH,jennite", "CSH,tobermorite", "SiO2(am)", "Calcite", "Al(OH)3(mic)", "Monosulfoaluminate", "Tricarboaluminate", "Monocarboaluminate", "Hemicarboaluminate", "Straetlingite", "Gypsum", "Ettringite", "Thaumasite" }; thedatabase.remove_gas_phases(); std::string co2_gas = R"plop( [ { "label": "CO2(g)", "composition": "CO2", "log_k": -1.468 } ] )plop"; thedatabase.add_gas_phases(co2_gas); specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); specmicp::index_t id_h2o = raw_data->get_id_component("H2O"); specmicp::index_t id_ho = raw_data->get_id_component("HO[-]"); specmicp::index_t id_co2 = raw_data->get_id_component("HCO3[-]"); specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = 1; specmicp::AdimensionalSystemSolverOptions options; options.solver_options.maxstep = 20.0; options.solver_options.max_iter = 100; options.solver_options.maxiter_maxstep = 100; options.solver_options.use_crashing = false; options.solver_options.use_scaling = true; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 100; options.solver_options.projection_min_variable = 1e-9; options.solver_options.fvectol = 1e-6; options.solver_options.steptol = 1e-10; options.system_options.non_ideality_tolerance = 1e-12; specmicp::Vector x; specmicp::scalar_t dh2co3 = mult*0.01; specmicp::uindex_t tot_iter = 0; specmicp::Timer timer_carboalu; timer_carboalu.start(); specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); - solver.initialise_variables(x, 0.5, -3.0); + solver.initialize_variables(x, 0.5, -3.0); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false); REQUIRE((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); std::cout << "First point number of iterations : " << perf.nb_iterations << " (Ref: 38)"<< std::endl; tot_iter += perf.nb_iterations; specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); for (int k=0; k<250; ++k) { constraints.total_concentrations(id_h2o) += dh2co3; constraints.total_concentrations(id_ho) -= dh2co3; constraints.total_concentrations(id_co2) += dh2co3; solver = specmicp::AdimensionalSystemSolver (raw_data, constraints, solution, options); //std::cout << solver.get_options().solver_options.factor_descent_condition << std::endl; specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false); //std::cout << x << std::endl; REQUIRE((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); //std::cout << perf.nb_iterations << std::endl; tot_iter += perf.nb_iterations; solution = solver.get_raw_solution(x); //specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); //std::cout << solution.main_variables(0) << std::endl; //std::cout << 14+solution.main_variables(1) << std::endl; } std::cout << "Number of iterations : " << tot_iter << " (Ref: 443)"<< std::endl; timer_carboalu.stop(); std::cout << "Elapsed time : " << timer_carboalu.elapsed_time() << "s (Ref: 0.022s)" << std::endl; } } diff --git a/tests/specmicp/adim/adimensional_system_conditions.cpp b/tests/specmicp/adim/adimensional_system_conditions.cpp index ba2596a..c97b29b 100644 --- a/tests/specmicp/adim/adimensional_system_conditions.cpp +++ b/tests/specmicp/adim/adimensional_system_conditions.cpp @@ -1,376 +1,376 @@ #include "catch.hpp" #include "specmicp_common/log.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp_database/database.hpp" #include TEST_CASE("Fixed activity", "[specmicp, MiCP, program, adimensional, solver, conditions]") { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Error; SECTION("Thermocarbo - saturated ") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); thedatabase.remove_gas_phases(); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); thedatabase.remove_half_cell_reactions({"H2O", "HO[-]", "HCO3[-]"}); specmicp::Formulation formulation; specmicp::scalar_t mult = 1.0; formulation.mass_solution = mult*0.156; formulation.amount_minerals = {{"C3S", mult*0.7}, {"C2S", mult*0.3}}; formulation.extra_components_to_keep = {"HCO3[-]", }; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = 1; constraints.water_equation = specmicp::WaterEquationType::SaturatedSystem; specmicp::index_t id_h2o = raw_data->water_index(); specmicp::index_t id_ho = raw_data->get_id_component("HO[-]"); specmicp::index_t id_co2 = raw_data->get_id_component("HCO3[-]"); specmicp::AdimensionalSystemSolverOptions options; options.solver_options.maxstep = 20.0; options.solver_options.maxiter_maxstep = 100; options.solver_options.use_crashing = false; options.solver_options.use_scaling = true; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 200; options.solver_options.condition_limit = -1; options.solver_options.set_tolerance(1e-10,1e-12); options.system_options.cutoff_total_concentration = 1e-16; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); specmicp::Vector x; - solver.initialise_variables(x, 0.8, 2); + solver.initialize_variables(x, 0.8, 2); specmicp::scalar_t dh2co3 = 0.1; for (int k=0; k<40; ++k) { specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor sol(solution, raw_data, options.units_set); //std::cout << solution.main_variables(0) << std::endl; std::cout << sol.pH() << std::endl; constraints.total_concentrations(id_h2o) += mult*dh2co3; constraints.total_concentrations(id_ho) -= mult*dh2co3; constraints.total_concentrations(id_co2) += mult*dh2co3; } } SECTION("Carbonate speciation - fixed activity") { specmicp::database::Database dbhandler(TEST_CEMDATA_PATH); std::vector to_keep = {"H2O", "H[+]", "HCO3[-]", "Na[+]"}; dbhandler.remove_gas_phases(); std::string co2_gas = R"plop( [ { "label": "CO2(g)", "composition": "CO2", "log_k": -1.468 } ] )plop"; dbhandler.remove_half_cell_reactions({"H2O", "H[+]", "HCO3[-]"}); dbhandler.add_gas_phases(co2_gas); dbhandler.keep_only_components(to_keep); specmicp::RawDatabasePtr thedatabase = dbhandler.get_database(); specmicp::index_t id_h = dbhandler.component_label_to_id("H[+]"); specmicp::index_t id_na = dbhandler.component_label_to_id("Na[+]"); specmicp::index_t id_hco3 = dbhandler.component_label_to_id("HCO3[-]"); specmicp::index_t id_co2 = dbhandler.aqueous_label_to_id("CO2"); specmicp::index_t id_co3 = dbhandler.aqueous_label_to_id("CO3[2-]"); specmicp::index_t id_naco3 = dbhandler.aqueous_label_to_id("NaCO3[-]"); specmicp::index_t id_nahco3 = dbhandler.aqueous_label_to_id("NaHCO3"); specmicp::index_t id_co2g = dbhandler.gas_label_to_id("CO2(g)"); specmicp::Vector total_concentrations(thedatabase->nb_component()); total_concentrations(0) = 55; total_concentrations(id_hco3) = 0.1; specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = id_na; constraints.enable_conservation_water(); specmicp::AdimensionalSystemSolverOptions options; options.system_options.non_ideality = true; options.units_set.length = specmicp::units::LengthUnit::decimeter; options.solver_options.max_iter = 50; options.solver_options.maxstep = 10; options.solver_options.maxiter_maxstep = 100; options.solver_options.use_crashing = false; options.solver_options.use_scaling = false; options.solver_options.non_monotone_linesearch = false; options.solver_options.penalization_factor = 0.8; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 10; options.solver_options.fvectol = 1e-6; specmicp::Vector x; constraints.add_fixed_activity_component(id_h, -4.0); specmicp::AdimensionalSystemSolver solver(thedatabase, constraints, options); - solver.initialise_variables(x, 0.8, -2); + solver.initialize_variables(x, 0.8, -2); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE((int) perf.return_code > 0 ); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set); std::cout << 4.0 << " \t" << sol.pH() << " \t" << sol.molality_component(id_hco3) << " \t" << sol.molality_aqueous(id_co2) << " \t" << sol.molality_aqueous(id_co3) << " \t" << sol.molality_aqueous(id_naco3) << " \t" << sol.molality_aqueous(id_nahco3) << "\t" << sol.fugacity_gas(id_co2g) << std::endl; for (double ph=4.5; ph<13; ph+=0.5) { constraints.fixed_activity_cs[0].log_value = -ph; solver = specmicp::AdimensionalSystemSolver( thedatabase, constraints, options); perf = solver.solve(x); REQUIRE((int) perf.return_code > 0 ); solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set); std::cout << ph << "\t" << sol.pH() << " \t" << sol.molality_component(id_hco3) << " \t" << sol.molality_aqueous(id_co2) << " \t" << sol.molality_aqueous(id_co3) << "\t" << sol.molality_aqueous(id_naco3) << "\t" << sol.molality_aqueous(id_nahco3) << "\t" << sol.fugacity_gas(id_co2g) << std::endl; } //std::cout << x << std::endl; } SECTION("Carbonate speciation - fixed molality") { specmicp::database::Database dbhandler(TEST_CEMDATA_PATH); std::vector to_keep = {"H2O", "H[+]", "HCO3[-]", "Na[+]", "Cl[-]"}; dbhandler.keep_only_components(to_keep); dbhandler.remove_gas_phases(); dbhandler.remove_half_cell_reactions({"H2O", "H[+]", "HCO3[-]"}); std::string co2_gas = R"plop( [ { "label": "CO2(g)", "composition": "CO2", "log_k": -1.468 } ] )plop"; dbhandler.add_gas_phases(co2_gas); specmicp::RawDatabasePtr thedatabase = dbhandler.get_database(); specmicp::index_t id_h = dbhandler.component_label_to_id("H[+]"); specmicp::index_t id_na = dbhandler.component_label_to_id("Na[+]"); specmicp::index_t id_cl = dbhandler.component_label_to_id("Cl[-]"); specmicp::index_t id_hco3 = dbhandler.component_label_to_id("HCO3[-]"); specmicp::index_t id_co2 = dbhandler.aqueous_label_to_id("CO2"); specmicp::index_t id_co3 = dbhandler.aqueous_label_to_id("CO3[2-]"); specmicp::index_t id_naco3 = dbhandler.aqueous_label_to_id("NaCO3[-]"); specmicp::index_t id_nahco3 = dbhandler.aqueous_label_to_id("NaHCO3"); specmicp::index_t id_co2g = dbhandler.gas_label_to_id("CO2(g)"); specmicp::Vector total_concentrations(thedatabase->nb_component()); total_concentrations(0) = 55; total_concentrations(id_hco3) = 0.1; specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.enable_conservation_water(); specmicp::AdimensionalSystemSolverOptions options; options.system_options.non_ideality = true; options.units_set.length = specmicp::units::LengthUnit::decimeter; options.solver_options.max_iter = 50; options.solver_options.maxstep = 10; options.solver_options.maxiter_maxstep = 100; options.solver_options.use_crashing = false; options.solver_options.use_scaling = true; options.solver_options.non_monotone_linesearch = true; options.solver_options.penalization_factor = 0.8; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 10; options.solver_options.fvectol = 1e-6; specmicp::Vector x; constraints.add_fixed_activity_component(id_h, -4.0); constraints.add_fixed_molality_component(id_na, -0.301); constraints.add_fixed_molality_component(id_cl, -0.301); //constraints.set_charge_keeper(id_h); specmicp::AdimensionalSystemSolver solver(thedatabase, constraints, options); - solver.initialise_variables(x, 0.8, -2); + solver.initialize_variables(x, 0.8, -2); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE((int) perf.return_code > 0 ); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set); REQUIRE(sol.log_molality_component(id_na) == Approx(-0.301)); REQUIRE(sol.log_molality_component(id_cl) == Approx(-0.301)); std::cout << 4.0 << " \t" << sol.pH() << " \t" << sol.molality_component(id_hco3) << " \t" << sol.molality_aqueous(id_co2) << " \t" << sol.molality_aqueous(id_co3) << " \t" << sol.molality_aqueous(id_naco3) << " \t" << sol.molality_aqueous(id_nahco3) << "\t" << sol.fugacity_gas(id_co2g) << std::endl; for (double ph=4.5; ph<13; ph+=0.5) { constraints.fixed_activity_cs[0].log_value = -ph; solver = specmicp::AdimensionalSystemSolver( thedatabase, constraints, options); perf = solver.solve(x); REQUIRE((int) perf.return_code > 0 ); solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set); std::cout << ph << "\t" << sol.pH() << " \t" << sol.molality_component(id_hco3) << " \t" << sol.molality_aqueous(id_co2) << " \t" << sol.molality_aqueous(id_co3) << "\t" << sol.molality_aqueous(id_naco3) << "\t" << sol.molality_aqueous(id_nahco3) << "\t" << sol.fugacity_gas(id_co2g) << std::endl; } //std::cout << x << std::endl; } SECTION("Carbonate speciation - fixed fugacity") { specmicp::database::Database dbhandler(TEST_CEMDATA_PATH); std::vector to_keep = {"H2O", "H[+]", "HCO3[-]"}; dbhandler.keep_only_components(to_keep); dbhandler.remove_gas_phases(); dbhandler.remove_half_cell_reactions({"H2O", "H[+]", "HCO3[-]"}); std::string co2_gas = R"plop( [ { "label": "CO2(g)", "composition": "CO2", "log_k": -1.468 } ] )plop"; dbhandler.add_gas_phases(co2_gas); specmicp::RawDatabasePtr thedatabase = dbhandler.get_database(); specmicp::index_t id_h = thedatabase->get_id_component("H[+]"); specmicp::index_t id_hco3 = thedatabase->get_id_component("HCO3[-]"); specmicp::index_t id_co2 = thedatabase->get_id_aqueous("CO2"); specmicp::index_t id_co3 = thedatabase->get_id_aqueous("CO3[2-]");; specmicp::index_t id_co2g = thedatabase->get_id_gas("CO2(g)"); CHECK(id_co2g == 0); specmicp::Vector total_concentrations(thedatabase->nb_component()); total_concentrations(0) = 55; total_concentrations(id_hco3) = 0.1; specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = id_h; constraints.enable_conservation_water(); specmicp::AdimensionalSystemSolverOptions options; options.system_options.non_ideality = true; options.units_set.length = specmicp::units::LengthUnit::decimeter; options.solver_options.max_iter = 50; options.solver_options.maxstep = 20; options.solver_options.maxiter_maxstep = 10; options.solver_options.use_crashing = false; options.solver_options.use_scaling = true; options.solver_options.enable_non_monotone_linesearch(); options.solver_options.penalization_factor = 0.8; options.solver_options.factor_descent_condition = -1; options.solver_options.fvectol = 1e-8; options.solver_options.condition_limit = -1; specmicp::Vector x; constraints.add_fixed_fugacity_gas(id_co2g, id_hco3, -5); specmicp::AdimensionalSystemSolver solver(thedatabase, constraints, options); - solver.initialise_variables(x, 0.8, -2.0); + solver.initialize_variables(x, 0.8, -2.0); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE((int) perf.return_code > 0 ); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set); std::cout << -5 << " \t" << sol.pH() << " \t" << sol.molality_component(id_hco3) << " \t" << sol.molality_aqueous(id_co2) << " \t" << sol.molality_aqueous(id_co3) << " \t" << sol.fugacity_gas(id_co2g) << std::endl; for (double lfug=-4.5; lfug<=-0.5; lfug+=0.25) { constraints.fixed_fugacity_cs[0].log_value = lfug; solver = specmicp::AdimensionalSystemSolver( thedatabase, constraints, options); perf = solver.solve(x); REQUIRE((int) perf.return_code > 0 ); solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor sol(solution, thedatabase, options.units_set); std::cout << lfug << "\t" << sol.pH() << " \t" << sol.molality_component(id_hco3) << " \t" << sol.molality_aqueous(id_co2) << " \t" << sol.molality_aqueous(id_co3) << "\t" << sol.fugacity_gas(id_co2g) << std::endl; } //std::cout << x << std::endl; } } diff --git a/tests/specmicp/adim/adimensional_system_extractor.cpp b/tests/specmicp/adim/adimensional_system_extractor.cpp index 38521bb..e0559d9 100644 --- a/tests/specmicp/adim/adimensional_system_extractor.cpp +++ b/tests/specmicp/adim/adimensional_system_extractor.cpp @@ -1,93 +1,93 @@ #include "catch.hpp" #include "specmicp_common/log.hpp" #include "specmicp/adimensional/adimensional_system.hpp" #include "specmicp_common/micpsolver/micpsolver.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp/adimensional/adimensional_system_solution_saver.hpp" #include "specmicp_database/database.hpp" #include static specmicp::RawDatabasePtr get_test_simple_database() { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, }); thedatabase.swap_components(swapping); std::vector to_keep = {"HO[-]", "Ca[2+]"}; thedatabase.keep_only_components(to_keep); thedatabase.remove_half_cell_reactions(std::vector({"H2O", "HO[-]",})) ; return thedatabase.get_database(); } using namespace specmicp; TEST_CASE("AdimSystemExtractor", "[extractor],[units]") { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Error; specmicp::RawDatabasePtr thedatabase = get_test_simple_database(); auto id_h2o = database::DataContainer::water_index(); auto id_oh = thedatabase->get_id_component("HO[-]"); auto id_ca = thedatabase->get_id_component("Ca[2+]"); auto id_ch = thedatabase->get_id_mineral("Portlandite"); SECTION("Automatic solver") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 0.03; total_concentration(id_oh) = 0.02; total_concentration(id_ca) = 0.01; specmicp::Vector x; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::AdimensionalSystemSolver solver(thedatabase, constraints); - solver.initialise_variables(x, 0.8, -2.0); + solver.initialize_variables(x, 0.8, -2.0); //x(solver.dof_surface()) = -HUGE_VAL; solver.get_options().units_set.length = specmicp::units::LengthUnit::centimeter; solver.get_options().solver_options.maxstep = 10.0; solver.get_options().solver_options.maxiter_maxstep = 100; solver.get_options().solver_options.use_crashing = false; solver.get_options().solver_options.use_scaling = true; solver.get_options().solver_options.disable_descent_direction(); solver.get_options().solver_options.factor_gradient_search_direction = 100; solver.solve(x); AdimensionalSystemSolution solution = solver.get_raw_solution(x); AdimensionalSystemSolutionExtractor extr(solution, thedatabase, solver.get_options().units_set); CHECK(extr.volume_fraction_water() == Approx(0.542049).epsilon(1e-4)); CHECK(extr.log_molality_component(id_oh) == Approx(-1.46214).epsilon(1e-4)); CHECK(extr.log_molality_component(id_ca) == Approx(-1.82093).epsilon(1e-4)); //CHECK(extr.free_surface_concentration() == -HUGE_VAL); CHECK(extr.volume_fraction_mineral(id_ch) == Approx(0.32966).epsilon(1e-4)); AdimensionalSystemSolutionSaver saver(thedatabase); saver.save_solution("test_solution.yaml", solution); AdimensionalSystemSolution solution2 = saver.parse_solution("test_solution.yaml"); REQUIRE(solution.main_variables(0) == Approx(solution2.main_variables(0))); REQUIRE(solution.main_variables(id_oh) == Approx(solution2.main_variables(id_oh))); REQUIRE(solution.main_variables(id_ca) == Approx(solution2.main_variables(id_ca))); REQUIRE(extr.volume_fraction_mineral(id_ch) == Approx(solution2.main_variables(extr.dof_mineral(id_ch)))); REQUIRE(solution.log_gamma.norm() == Approx(solution2.log_gamma.norm())); } } diff --git a/tests/specmicp/adim/adimensional_system_oxydoreduc.cpp b/tests/specmicp/adim/adimensional_system_oxydoreduc.cpp index a4bbff1..ebd898d 100644 --- a/tests/specmicp/adim/adimensional_system_oxydoreduc.cpp +++ b/tests/specmicp/adim/adimensional_system_oxydoreduc.cpp @@ -1,80 +1,80 @@ #include "catch.hpp" #include "specmicp_database/database.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include #include "specmicp_common/log.hpp" specmicp::RawDatabasePtr get_test_s_oxydoreduc_database() { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, }); thedatabase.swap_components(swapping); std::vector to_keep = {"HO[-]", "Ca[2+]"}; thedatabase.keep_only_components(to_keep); return thedatabase.get_database(); } using namespace specmicp; TEST_CASE("Oxydoreduction solver", "[Adimensional][OxydoReduc],[Solver]") { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Error; specmicp::RawDatabasePtr thedatabase = get_test_s_oxydoreduc_database(); auto id_h2o = database::DataContainer::water_index(); auto id_oh = thedatabase->get_id_component("HO[-]"); auto id_ca = thedatabase->get_id_component("Ca[2+]"); auto id_ch = thedatabase->get_id_mineral("Portlandite"); SECTION("Automatic solver") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 0.03; total_concentration(id_oh) = 0.02; total_concentration(id_ca) = 0.01; specmicp::Vector x; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::AdimensionalSystemSolver solver(thedatabase, constraints); - solver.initialise_variables(x, 0.8, -2.0); + solver.initialize_variables(x, 0.8, -2.0); //x(solver.dof_surface()) = -HUGE_VAL; solver.get_options().units_set.length = specmicp::units::LengthUnit::centimeter; solver.get_options().solver_options.maxstep = 10.0; solver.get_options().solver_options.maxiter_maxstep = 100; solver.get_options().solver_options.use_crashing = false; solver.get_options().solver_options.use_scaling = true; solver.get_options().solver_options.disable_descent_direction(); solver.get_options().solver_options.factor_gradient_search_direction = 100; solver.get_options().solver_options.set_tolerance(1e-8, 1e-14); solver.get_options().system_options.scaling_electron = 1e-8; solver.solve(x); AdimensionalSystemSolution solution = solver.get_raw_solution(x); AdimensionalSystemSolutionExtractor extr(solution, thedatabase, solver.get_options().units_set); CHECK(extr.volume_fraction_water() == Approx(0.542049).epsilon(1e-3)); CHECK(extr.log_molality_component(id_oh) == Approx(-1.46214).epsilon(1e-3)); CHECK(extr.log_molality_component(id_ca) == Approx(-1.82093).epsilon(1e-3)); //CHECK(extr.free_surface_concentration() == -HUGE_VAL); CHECK(extr.volume_fraction_mineral(id_ch) == Approx(0.32966).epsilon(1e-3)); CHECK(extr.pE() == Approx(1.32169).epsilon(1e-3)); CHECK(extr.Eh() == Approx(0.0775594).epsilon(1e-3)); } } diff --git a/tests/specmicp/adim/adimensional_system_pcfm.cpp b/tests/specmicp/adim/adimensional_system_pcfm.cpp index 238ca59..58f5c0f 100644 --- a/tests/specmicp/adim/adimensional_system_pcfm.cpp +++ b/tests/specmicp/adim/adimensional_system_pcfm.cpp @@ -1,61 +1,61 @@ #include "catch.hpp" #include "specmicp/adimensional/adimensional_system.hpp" #include "specmicp/adimensional/adimensional_system_pcfm.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp_database/database.hpp" specmicp::RawDatabasePtr get_pcfm_test_database() { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, }); thedatabase.swap_components(swapping); std::vector to_keep = {"H2O", "HO[-]", "Ca[2+]"}; thedatabase.keep_only_components(to_keep); thedatabase.remove_solid_phases(); return thedatabase.get_database(); } // Test the Positive continuous fraction method TEST_CASE("Positive continuous fraction method", "[specmicp, MiCP, program, adimensional, PCFM]") { specmicp::RawDatabasePtr the_database = get_pcfm_test_database(); auto id_h2o = the_database->water_index(); auto id_oh = the_database->get_id_component("HO[-]"); auto id_ca = the_database->get_id_component("Ca[2+]"); SECTION("PCFM") { specmicp::Vector total_concentration = specmicp::Vector::Zero(the_database->nb_component()); total_concentration(id_h2o) = 55.5; total_concentration(id_oh) = 2e-3; total_concentration(id_ca) = 1e-3; specmicp::AdimensionalSystemConstraints constraints(total_concentration); constraints.disable_conservation_water(); constraints.enable_surface_model(1.23456); //constraints.total_concentrations(2) = 0; std::shared_ptr ptrsystem = std::make_shared(the_database, constraints); specmicp::AdimensionalSystemSolver solver(the_database, constraints); specmicp::AdimensionalSystemPCFM pcfm_solver(ptrsystem); specmicp::Vector x; - solver.initialise_variables(x, 1.0, -3.0, 0.0); + solver.initialize_variables(x, 1.0, -3.0, 0.0); specmicp::PCFMReturnCode retcode = pcfm_solver.solve(x); REQUIRE(retcode == specmicp::PCFMReturnCode::Success); } } diff --git a/tests/specmicp/adim/adimensional_system_problem_solver.cpp b/tests/specmicp/adim/adimensional_system_problem_solver.cpp index ec123ac..063e929 100644 --- a/tests/specmicp/adim/adimensional_system_problem_solver.cpp +++ b/tests/specmicp/adim/adimensional_system_problem_solver.cpp @@ -1,368 +1,369 @@ #include "catch.hpp" #include "specmicp_common/log.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "specmicp/problem_solver/reactant_box.hpp" +#include "specmicp/problem_solver/smart_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp_database/database.hpp" #include TEST_CASE("Solving adimensinal system with problem setting", "[specmicp, MiCP, program, adimensional, solver]") { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; SECTION("Solving problem with dissolver - simpler problem") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::Formulation formulation; formulation.mass_solution = 0.8; formulation.amount_minerals = {{"Portlandite", 10.0}}; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); total_concentrations *= 1e3; specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = raw_data->get_id_component("HO[-]"); specmicp::AdimensionalSystemSolver solver(raw_data, constraints); solver.get_options().solver_options.maxstep = 50.0; solver.get_options().solver_options.maxiter_maxstep = 100; solver.get_options().solver_options.use_crashing = false; solver.get_options().solver_options.use_scaling = true; solver.get_options().solver_options.factor_descent_condition = 1e-10; solver.get_options().solver_options.factor_gradient_search_direction = 200; specmicp::Vector x; - solver.initialise_variables(x, 0.8, -2); + solver.initialize_variables(x, 0.8, -2); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(solution, raw_data, solver.get_options().units_set); CHECK(extr.volume_fraction_water() == Approx(0.802367).epsilon(1e-4)); CHECK(extr.volume_fraction_mineral(raw_data->get_id_mineral("Portlandite")) == Approx(0.32947).epsilon(1e-2)); } SECTION("Solving problem with dissolver") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions({"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::Formulation formulation; formulation.mass_solution = 0.8; formulation.amount_minerals = {{"C3S", 0.7}, {"C2S", 0.3}}; formulation.extra_components_to_keep = {"HCO3[-]", }; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); total_concentrations *= 1e3; total_concentrations(2) = 500; specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = raw_data->get_id_component("HO[-]"); specmicp::AdimensionalSystemSolver solver(raw_data, constraints); solver.get_options().solver_options.maxstep = 20.0; solver.get_options().solver_options.maxiter_maxstep = 100; solver.get_options().solver_options.use_crashing = false; solver.get_options().solver_options.use_scaling = true; solver.get_options().solver_options.factor_descent_condition = 1e-10; solver.get_options().solver_options.factor_gradient_search_direction = 200; specmicp::Vector x(raw_data->nb_component()+raw_data->nb_mineral()); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, true); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); } SECTION("Solving problem with reactant box") { std::cout << "\n Solving with reactant box \n ------------------------ \n"; std::cout << "== with SI units == \n"; specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"}, {"HCO3[-]", "CO2"} }); thedatabase.swap_components(swapping); //thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(500, "L"); reactant_box.add_solid_phase("C3S", 350, "L"); reactant_box.add_solid_phase("C2S", 50, "L"); reactant_box.add_aqueous_species("CO2", 32.72, "mol/kg"); reactant_box.set_charge_keeper("HO[-]"); specmicp::AdimensionalSystemConstraints constraints = reactant_box.get_constraints(true); specmicp::Vector total_concentrations = constraints.total_concentrations; specmicp::AdimensionalSystemSolver solver(raw_data, constraints); specmicp::Vector x; - solver.initialise_variables(x, 0.5, -6); + solver.initialize_variables(x, 0.5, -6); solver.get_options().solver_options.set_tolerance(1e-8, 1e-10); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); std::cout << "perf nb it : " << perf.nb_iterations << std::endl; REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(solution,raw_data, specmicp::units::SI_units); std::cout << "== with non SI units : mmol/dm == \n"; specmicp::units::UnitsSet dmmol; dmmol.length = specmicp::units::LengthUnit::decimeter; dmmol.quantity = specmicp::units::QuantityUnit::millimoles; reactant_box.set_units(dmmol); specmicp::AdimensionalSystemConstraints constraints2 = reactant_box.get_constraints(false); specmicp::Vector total_concentrations2 = constraints2.total_concentrations; specmicp::AdimensionalSystemSolver solver2(raw_data, constraints2); solver2.get_options().units_set = dmmol; solver2.get_options().solver_options.set_tolerance(1e-8, 1e-10); specmicp::Vector x2; - solver2.initialise_variables(x2, 0.5, -6); + solver2.initialize_variables(x2, 0.5, -6); perf = solver2.solve(x2); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); std::cout << "perf nb it : " << perf.nb_iterations << std::endl; specmicp::AdimensionalSystemSolution solution2 = solver2.get_raw_solution(x); CHECK(solution.main_variables(0) == Approx(solution2.main_variables(0)).epsilon(1e-6)); CHECK(solution.main_variables(2) == Approx(solution2.main_variables(2)).epsilon(1e-6)); CHECK(solution.main_variables(5) == Approx(solution2.main_variables(5)).epsilon(1e-6)); CHECK(solution.ionic_strength == Approx(solution2.ionic_strength).epsilon(1e-6)); specmicp::AdimensionalSystemSolutionExtractor extr2(solution2, raw_data, dmmol); for (auto idc: raw_data->range_component()) { if (idc == specmicp::database::electron_id) continue; // std::cout << "Id " << idc << " - " << raw_data->get_label_component(idc) << std::endl; // std::cout << total_concentrations(idc) << " - " // << total_concentrations2(idc) << " - " // << extr.total_concentration(idc) << " ~ " // << extr.volume_fraction_gas_phase() * extr.total_gaseous_concentration(idc) << " # " // << extr.fugacity_gas(0) << " - " // << extr2.total_concentration(idc) // << std::endl; CHECK(extr.total_concentration(idc) == Approx(total_concentrations(idc)).epsilon(1e-8)); CHECK(extr2.total_concentration(idc) == Approx(total_concentrations2(idc)).epsilon(1e-8)); CHECK(extr.total_concentration(idc) == Approx(extr2.total_concentration(idc)).epsilon(1e-8)); } std::cout << "== with non SI units : mmol/cm == \n"; specmicp::units::UnitsSet cmmol; cmmol.length = specmicp::units::LengthUnit::centimeter; cmmol.quantity = specmicp::units::QuantityUnit::millimoles; reactant_box.set_units(cmmol); specmicp::AdimensionalSystemConstraints constraints3 = reactant_box.get_constraints(false); specmicp::Vector total_concentrations3 = constraints3.total_concentrations; specmicp::AdimensionalSystemSolver solver3(raw_data, constraints3); solver3.get_options().units_set = cmmol; solver3.get_options().solver_options.set_tolerance(1e-8, 1e-12); specmicp::Vector x3; - solver3.initialise_variables(x3, 0.5, -6); + solver3.initialize_variables(x3, 0.5, -6); perf = solver3.solve(x3); for (auto idc: raw_data->range_component()) { if (idc == specmicp::database::electron_id) continue; CHECK(total_concentrations(idc) == Approx(1e3*total_concentrations3(idc)).epsilon(1e-8)); } REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); std::cout << "perf nb it : " << perf.nb_iterations << std::endl; specmicp::AdimensionalSystemSolution solution3 = solver3.get_raw_solution(x); CHECK(solution.main_variables(0) == Approx(solution3.main_variables(0)).epsilon(1e-6)); CHECK(solution.main_variables(2) == Approx(solution3.main_variables(2)).epsilon(1e-6)); CHECK(solution.main_variables(5) == Approx(solution3.main_variables(5)).epsilon(1e-6)); CHECK(solution.ionic_strength == Approx(solution3.ionic_strength).epsilon(1e-6)); specmicp::AdimensionalSystemSolutionExtractor extr3(solution3, raw_data, cmmol); for (auto idc: raw_data->range_component()) { CHECK(extr3.total_concentration(idc) == Approx(total_concentrations3(idc)).epsilon(1e-8)); } } SECTION("Reactant box fixed activity") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); //thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(500, "L"); reactant_box.add_solid_phase("C3S", 350, "L"); reactant_box.add_solid_phase("C2S", 50, "L"); reactant_box.add_aqueous_species("CO2", 15, "mol/kg"); reactant_box.set_charge_keeper("Cl[-]"); reactant_box.add_fixed_activity_component("H[+]", 1e-8); auto constraints = reactant_box.get_constraints(true); specmicp::AdimensionalSystemSolver solver(raw_data, constraints); specmicp::Vector x; - solver.initialise_variables(x, 0.8, -6); + solver.initialize_variables(x, 0.8, -6); solver.get_options().solver_options.set_tolerance(1e-8, 1e-10); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units); CHECK(extr.pH() == Approx(8)); } SECTION("Reactant box fixed molality") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); //thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(500, "L"); reactant_box.add_solid_phase("C3S", 350, "L"); reactant_box.add_solid_phase("C2S", 50, "L"); reactant_box.add_aqueous_species("CO2", 15, "mol/kg"); reactant_box.set_charge_keeper("Cl[-]"); reactant_box.add_fixed_molality_component("H[+]", 1e-8); auto constraints = reactant_box.get_constraints(true); specmicp::AdimensionalSystemSolver solver(raw_data, constraints); specmicp::Vector x; - solver.initialise_variables(x, 0.8, -6); + solver.initialize_variables(x, 0.8, -6); solver.get_options().solver_options.set_tolerance(1e-8, 1e-10); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units); CHECK(extr.molality_component(raw_data->get_id_component("H[+]")) == Approx(1e-8)); } SECTION("Reactant box fixed fugacity") { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); //thedatabase.remove_gas_phases(); thedatabase.remove_half_cell_reactions();//{"H2O", "HO[-]", "HCO3[-]"}); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); specmicp::ReactantBox reactant_box(raw_data, specmicp::units::SI_units); reactant_box.set_solution(500, "L"); reactant_box.add_solid_phase("C3S", 350, "L"); reactant_box.add_solid_phase("C2S", 50, "L"); reactant_box.add_fixed_fugacity_gas("CO2(g)", "HCO3[-]", 1e-3); reactant_box.set_charge_keeper("H[+]"); auto constraints = reactant_box.get_constraints(true); specmicp::AdimensionalSystemSolver solver(raw_data, constraints); specmicp::Vector x; - solver.initialise_variables(x, 0.8, -4); + solver.initialize_variables(x, 0.8, -4); solver.get_options().solver_options.set_tolerance(1e-6, 1e-8); specmicp::micpsolver::MiCPPerformance perf = solver.solve(x); REQUIRE(perf.return_code >= specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet); specmicp::AdimensionalSystemSolution sol = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(sol, raw_data, specmicp::units::SI_units); CHECK(extr.fugacity_gas(raw_data->get_id_gas("CO2(g)")) == Approx(1e-3)); } } diff --git a/tests/specmicp/adim/adimensional_system_solver.cpp b/tests/specmicp/adim/adimensional_system_solver.cpp index 6c5a169..800b97c 100644 --- a/tests/specmicp/adim/adimensional_system_solver.cpp +++ b/tests/specmicp/adim/adimensional_system_solver.cpp @@ -1,221 +1,221 @@ #include "catch.hpp" #include "specmicp_common/log.hpp" #include "specmicp/adimensional/adimensional_system.hpp" #include "specmicp_common/micpsolver/micpsolver.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp/adimensional/adimensional_system_solution_saver.hpp" #include "specmicp_database/database.hpp" #include specmicp::RawDatabasePtr get_test_simple_database() { specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); std::map swapping ({ {"H[+]","HO[-]"}, }); thedatabase.swap_components(swapping); std::vector to_keep = {"HO[-]", "Ca[2+]"}; thedatabase.keep_only_components(to_keep); thedatabase.remove_half_cell_reactions(std::vector({"H2O", "HO[-]",})) ; return thedatabase.get_database(); } using namespace specmicp; TEST_CASE("Solving adimensional system", "[specmicp, MiCP, program, adimensional, solver]") { specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Error; specmicp::RawDatabasePtr thedatabase = get_test_simple_database(); auto id_h2o = database::DataContainer::water_index(); auto id_oh = thedatabase->get_id_component("HO[-]"); auto id_ca = thedatabase->get_id_component("Ca[2+]"); auto id_ch = thedatabase->get_id_mineral("Portlandite"); SECTION("Solving simple case") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 3.0e4; total_concentration(id_oh) = 2.0e4; total_concentration(id_ca) = 1.0e4; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::micpsolver::MiCPSolverOptions options; options.maxstep = 100; options.maxiter_maxstep = 100; options.disable_crashing(); options.enable_scaling(); options.disable_descent_direction(); std::shared_ptr system= std::make_shared(thedatabase, constraints); specmicp::micpsolver::MiCPSolver solver(system); solver.set_options(options); Vector x = Vector::Zero(system->total_variables()); x(system->ideq_w()) = 0.8; x(system->ideq_paq(id_oh)) = -2.0; x(system->ideq_paq(id_ca)) = -2.3; system->compute_log_gamma(x); system->set_secondary_concentration(x); solver.solve(x); CHECK(x(system->ideq_w()) == Approx(0.542049).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.462142).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.820933).epsilon(1e-4)); CHECK(x(system->ideq_min(id_ch)) == Approx(0.32966).epsilon(1e-4)); } SECTION("Solving simple case - volume in cm^3") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 0.03; total_concentration(id_oh) = 0.02; total_concentration(id_ca) = 0.01; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::micpsolver::MiCPSolverOptions options; options.maxstep = 10; options.maxiter_maxstep = 100; options.disable_crashing(); options.enable_scaling(); options.disable_descent_direction(); std::shared_ptr system= std::make_shared(thedatabase, constraints); units::UnitsSet unit_set; unit_set.length = units::LengthUnit::centimeter; system->set_units(unit_set); specmicp::micpsolver::MiCPSolver solver(system); solver.set_options(options); Vector x = Vector::Zero(system->total_variables()); x(system->ideq_w()) = 0.8; x(system->ideq_paq(id_oh)) = -2.0; x(system->ideq_paq(id_ca)) = -2.3; system->compute_log_gamma(x); system->set_secondary_concentration(x); solver.solve(x); CHECK(x(system->ideq_w()) == Approx( 0.54205).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.46214).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.82093).epsilon(1e-4)); CHECK(x(system->ideq_min(id_ch)) == Approx( 0.32966).epsilon(1e-4)); } SECTION("Solving simple case - volume in cm^3, mmol") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 1e3*0.03; total_concentration(id_oh) = 1e3*0.02; total_concentration(id_ca) = 1e3*0.01; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::micpsolver::MiCPSolverOptions options; options.maxstep = 10; options.maxiter_maxstep = 100; options.disable_crashing(); options.enable_scaling(); options.disable_descent_direction(); std::shared_ptr system= std::make_shared(thedatabase, constraints); units::UnitsSet unit_set; unit_set.length = units::LengthUnit::centimeter; unit_set.quantity = units::QuantityUnit::millimoles; system->set_units(unit_set); specmicp::micpsolver::MiCPSolver solver(system); solver.set_options(options); Vector x = Vector::Zero(system->total_variables()); x(system->ideq_w()) = 0.8; x(system->ideq_paq(id_oh)) = -2.0; x(system->ideq_paq(id_ca)) = -2.3; system->compute_log_gamma(x); system->set_secondary_concentration(x); solver.solve(x); CHECK(x(system->ideq_w()) == Approx( 0.54205).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_oh)) == Approx(-1.46214).epsilon(1e-4)); CHECK(x(system->ideq_paq(id_ca)) == Approx(-1.82093).epsilon(1e-4)); CHECK(x(system->ideq_min(id_ch)) == Approx(0.32966).epsilon(1e-4)); } SECTION("Automatic solver") { Vector total_concentration = Vector::Zero(thedatabase->nb_component()); total_concentration(id_h2o) = 0.03; total_concentration(id_oh) = 0.02; total_concentration(id_ca) = 0.01; specmicp::Vector x; specmicp::AdimensionalSystemConstraints constraints(total_concentration); specmicp::AdimensionalSystemSolver solver(thedatabase, constraints); - solver.initialise_variables(x, 0.8, -2.0); + solver.initialize_variables(x, 0.8, -2.0); //x(solver.dof_surface()) = -HUGE_VAL; solver.get_options().units_set.length = specmicp::units::LengthUnit::centimeter; solver.get_options().solver_options.maxstep = 10.0; solver.get_options().solver_options.maxiter_maxstep = 100; solver.get_options().solver_options.use_crashing = false; solver.get_options().solver_options.use_scaling = true; solver.get_options().solver_options.disable_descent_direction(); solver.get_options().solver_options.factor_gradient_search_direction = 100; solver.solve(x); AdimensionalSystemSolution solution = solver.get_raw_solution(x); AdimensionalSystemSolutionExtractor extr(solution, thedatabase, solver.get_options().units_set); CHECK(extr.volume_fraction_water() == Approx(0.542049).epsilon(1e-4)); CHECK(extr.log_molality_component(id_oh) == Approx(-1.46214).epsilon(1e-4)); CHECK(extr.log_molality_component(id_ca) == Approx(-1.82093).epsilon(1e-4)); //CHECK(extr.free_surface_concentration() == -HUGE_VAL); CHECK(extr.volume_fraction_mineral(id_ch) == Approx(0.32966).epsilon(1e-4)); AdimensionalSystemSolutionSaver saver(thedatabase); saver.save_solution("test_solution.yaml", solution); AdimensionalSystemSolution solution2 = saver.parse_solution("test_solution.yaml"); REQUIRE(solution.main_variables(0) == Approx(solution2.main_variables(0))); REQUIRE(solution.main_variables(id_oh) == Approx(solution2.main_variables(id_oh))); REQUIRE(solution.main_variables(id_ca) == Approx(solution2.main_variables(id_ca))); REQUIRE(extr.volume_fraction_mineral(id_ch) == Approx(solution2.main_variables(extr.dof_mineral(id_ch)))); REQUIRE(solution.log_gamma.norm() == Approx(solution2.log_gamma.norm())); } } diff --git a/tests/specmicp/adim/adimensional_system_thermocarbo.cpp b/tests/specmicp/adim/adimensional_system_thermocarbo.cpp index 4f6d6fd..2cacab0 100644 --- a/tests/specmicp/adim/adimensional_system_thermocarbo.cpp +++ b/tests/specmicp/adim/adimensional_system_thermocarbo.cpp @@ -1,156 +1,156 @@ #include "catch.hpp" #include "specmicp_common/log.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "specmicp_common/timer.hpp" #include "specmicp_database/database.hpp" #include TEST_CASE("thermocarbo - using adimensional system ", "[Adimensional, Thermocarbo]") { std::cerr.flush(); std::cout.flush(); specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; SECTION("Thermocarbo") { std::cout << "\n-------------------\n Thermocarbo\n -------------------" << std::endl; specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); thedatabase.remove_gas_phases(); specmicp::Formulation formulation; specmicp::scalar_t mult = 2.0; formulation.mass_solution = mult*0.156; formulation.amount_minerals = {{"C3S", mult*0.7}, {"C2S", mult*0.3}}; formulation.extra_components_to_keep = {"HCO3[-]", }; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = raw_data->get_id_component("HO[-]"); specmicp::index_t id_h2o = raw_data->water_index(); specmicp::index_t id_ho = raw_data->get_id_component("HO[-]"); specmicp::index_t id_co2 = raw_data->get_id_component("HCO3[-]"); specmicp::AdimensionalSystemSolverOptions options; options.solver_options.maxstep = 20.0; options.solver_options.maxiter_maxstep = 100; options.solver_options.disable_descent_direction(); options.solver_options.set_tolerance(1e-8, 1e-12); //options.solver_options.disable_non_monotone_linesearch(); options.solver_options.disable_condition_check(); options.solver_options.disable_crashing(); options.solver_options.enable_scaling(); options.units_set.length = specmicp::units::LengthUnit::decimeter; specmicp::Timer tot_timer; tot_timer.start(); specmicp::Vector x; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); - solver.initialise_variables(x, 0.8, -4.0); + solver.initialize_variables(x, 0.8, -4.0); specmicp::scalar_t dh2co3 = 0.1; specmicp::index_t nb_iter {0}; constexpr specmicp::index_t nb_step = 40; specmicp::Vector solution(nb_step); for (int k=0; k(perf.return_code) >= static_cast(specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet)); nb_iter += perf.nb_iterations; constraints.total_concentrations(id_h2o) += mult*dh2co3; constraints.total_concentrations(id_ho) -= mult*dh2co3; constraints.total_concentrations(id_co2) += mult*dh2co3; specmicp::AdimensionalSystemSolution adim_solution = solver.get_raw_solution(x); solution(k) = 14+adim_solution.main_variables(id_ho); specmicp::AdimensionalSystemSolutionExtractor extr(adim_solution, raw_data, options.units_set); //std::cout << "Saturation water : " << extr.volume_fraction_water() << " - " // << extr.porosity() << " - " // << extr.saturation_water() << std::endl; } tot_timer.stop(); std::cout << "Total time : " << tot_timer.elapsed_time() << "s (Ref 0.01s)" << std::endl; std::cout << "Nb iter : " << nb_iter << " (Ref 560)" << std::endl; specmicp::Vector reference_solution(nb_step); reference_solution << 12.53785, 12.53785, 12.53785, 12.53785, 12.53785, 12.53785, 12.53785, 12.53785, 12.53785, 12.53785, 12.53785, 12.15348, 12.15348, 12.15348, 12.15348, 12.15348, 12.15348, 12.15348, 12.15348, 9.8575, 9.8575, 9.8575, 9.8575, 9.8575, 9.8575, 9.8575, 9.8575, 8.97554, 5.35498, 5.16745, 5.05979, 4.98442, 4.92662, 4.87987, 4.84073, 4.80714, 4.77778, 4.75174, 4.72839, 4.70725; for (int k=0;k using namespace specmicp; scalar_t water_pressure(scalar_t sat) { const scalar_t tmp = std::pow(sat, -2.1684) -1; const scalar_t exp = 1.0-1.0/2.1684; scalar_t ln_hr = - 3.67 * std::pow(tmp, exp); const scalar_t p = 31.7e2*std::exp(ln_hr); return p; // in Pa } TEST_CASE("thermocarbo - with water pressure", "[Adimensional],[Thermocarbo],[Water pressure]") { std::cerr.flush(); std::cout.flush(); specmicp::logger::ErrFile::stream() = &std::cerr; specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; SECTION("Thermocarbo") { std::cout << "\n-------------------\n Thermocarbo - with water pressure\n -------------------" << std::endl; specmicp::database::Database thedatabase(TEST_CEMDATA_PATH); specmicp::RawDatabasePtr raw_data = thedatabase.get_database(); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"} }); thedatabase.swap_components(swapping); thedatabase.remove_gas_phases(); specmicp::Formulation formulation; specmicp::scalar_t mult = 3.0; formulation.mass_solution = mult*0.156; formulation.amount_minerals = {{"C3S", mult*0.7}, {"C2S", mult*0.3}}; formulation.extra_components_to_keep = {"HCO3[-]", }; specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation); specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.charge_keeper = raw_data->get_id_component("HO[-]"); constraints.set_water_partial_pressure_model((water_partial_pressure_f) water_pressure); specmicp::index_t id_h2o = raw_data->water_index(); specmicp::index_t id_ho = raw_data->get_id_component("HO[-]"); specmicp::index_t id_co2 = raw_data->get_id_component("HCO3[-]"); specmicp::AdimensionalSystemSolverOptions options; options.solver_options.maxstep = 20.0; options.solver_options.maxiter_maxstep = 100; options.solver_options.disable_descent_direction(); options.solver_options.set_tolerance(1e-10, 1e-12); //options.solver_options.disable_non_monotone_linesearch(); options.solver_options.disable_condition_check(); options.solver_options.disable_crashing(); options.solver_options.enable_scaling(); options.units_set.length = specmicp::units::LengthUnit::decimeter; specmicp::Timer tot_timer; tot_timer.start(); specmicp::Vector x; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); - solver.initialise_variables(x, 0.8, -4.0); + solver.initialize_variables(x, 0.8, -4.0); specmicp::scalar_t dh2co3 = 0.1; specmicp::index_t nb_iter {0}; constexpr specmicp::index_t nb_step = 20; for (int k=0; k(perf.return_code) >= static_cast(specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet)); nb_iter += perf.nb_iterations; specmicp::AdimensionalSystemSolution adim_solution = solver.get_raw_solution(x); specmicp::AdimensionalSystemSolutionExtractor extr(adim_solution, raw_data, options.units_set); //std::cout << "Saturation water : " << extr.volume_fraction_water() << " / " << x(0) << " = " << extr.volume_fraction_water() / x(0) // << " - " // << extr.porosity() << " - " // << extr.saturation_water() << std::endl; const scalar_t tot_conc_w = extr.total_solid_concentration(0) + extr.mass_concentration_water()*extr.total_aqueous_concentration(0); const scalar_t tot_conc_g = extr.volume_fraction_gas_phase()*1e-3*( water_pressure(extr.saturation_water())/ (constants::gas_constant*(273.16+25.0)) ); const scalar_t tot_conc_w_pg = tot_conc_w + tot_conc_g; //std::cout << "Vol frac gas : " << extr.volume_fraction_gas_phase() << std::endl; //std::cout << "hop : " << extr.total_solid_concentration(0) << " + " // << extr.mass_concentration_water()*extr.total_aqueous_concentration(0) << " + " // << tot_conc_g << std::endl; //std::cout << "tot_conc : " << tot_conc_w0 << " - " << tot_conc_w << " = " << tot_conc_w0 - tot_conc_w << std::endl; //std::cout << "tot_conc : " << tot_conc_w0 << " - " << tot_conc_w_pg << " = " << tot_conc_w0 - tot_conc_w_pg << std::endl; CHECK(tot_conc_w0 == Approx(tot_conc_w_pg).epsilon(1e-10)); //std::cout << "diff tot " << tot_conc_w - tot_conc_w_pg << std::endl; //std::cout << "p_v(" << extr.saturation_water() << ") = " << water_pressure(extr.saturation_water()) << " Pa"<< std::endl; constraints.total_concentrations(id_h2o) += mult*dh2co3; constraints.total_concentrations(id_ho) -= mult*dh2co3; constraints.total_concentrations(id_co2) += mult*dh2co3; } tot_timer.stop(); std::cout << "Total time : " << tot_timer.elapsed_time() << "s (Ref 0.007s)" << std::endl; std::cout << "Nb iter : " << nb_iter << " (Ref 360)" << std::endl; } }