diff --git a/examples/reactmicp/saturated_react/momas_benchmark.cpp b/examples/reactmicp/saturated_react/momas_benchmark.cpp index 2f81f52..64f4dee 100644 --- a/examples/reactmicp/saturated_react/momas_benchmark.cpp +++ b/examples/reactmicp/saturated_react/momas_benchmark.cpp @@ -1,598 +1,604 @@ /* ============================================================================= 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 IS_ADVECTIVE false #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.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; + Vector sc(1); + sc(0) = 0.25*1.0; + constraints_a.surface_model.total_ssite_concentrations = sc; 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; + Vector sc(1); + sc(0) = 0.5*10.0; + constraints_b.surface_model.total_ssite_concentrations = sc; 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; + Vector sc(1); + sc(0) = 0.0; + constraints_bc.surface_model.total_ssite_concentrations = sc; 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.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.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.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; + //options.system_options.restart_concentration = -2; + //options.use_pcfm = true; return options; } mesh::Mesh1DPtr get_mesh(scalar_t dx) { mesh::Uniform1DMeshGeometry geom; scalar_t tot_length = 2.1+dx; geom.nb_nodes = tot_length/dx + 2; geom.dx = dx; geom.section = 1.0; return mesh::uniform_mesh1d(geom); } 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_x1 << "\t" << variables->aqueous_concentration(node, 2, variables->displacement()); + output_x2 << "\t" << variables->aqueous_concentration(node, 3, variables->displacement()); + output_x3 << "\t" << variables->aqueous_concentration(node, 4, variables->displacement()); + output_x4 << "\t" << variables->aqueous_concentration(node, 5, variables->displacement()); + output_sx1 << "\t" << variables->solid_concentration(node, 2, variables->displacement()); + output_sx2 << "\t" << variables->solid_concentration(node, 3, variables->displacement()); + output_sx3 << "\t" << variables->solid_concentration(node, 4, variables->displacement()); + output_sx4 << "\t" << variables->solid_concentration(node, 5, 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_s << "\t" << pow10(variables->equilibrium_solution(node).main_variables(6)); } 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_x1 << "\t" << variables->aqueous_concentration(node, 2, variables->displacement()); + output_x2 << "\t" << variables->aqueous_concentration(node, 3, variables->displacement()); + output_x3 << "\t" << variables->aqueous_concentration(node, 4, variables->displacement()); + output_x4 << "\t" << variables->aqueous_concentration(node, 5, variables->displacement()); + output_sx1 << "\t" << variables->solid_concentration(node, 2, variables->displacement()); + output_sx2 << "\t" << variables->solid_concentration(node, 3, variables->displacement()); + output_sx3 << "\t" << variables->solid_concentration(node, 4, variables->displacement()); + output_sx4 << "\t" << variables->solid_concentration(node, 5, 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_s << "\t" << pow10(variables->equilibrium_solution(node).main_variables(6)); } 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/specmicp/adimensional/momas_thermo.cpp b/examples/specmicp/adimensional/momas_thermo.cpp index 012d25b..239ef88 100644 --- a/examples/specmicp/adimensional/momas_thermo.cpp +++ b/examples/specmicp/adimensional/momas_thermo.cpp @@ -1,284 +1,288 @@ /* ============================================================================= 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; + std::cout << raw_data->nb_component() << std::endl; + total_concentrations << 1, 0.0, 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::Vector sc(1); + sc(0) = 0.5; + constraints().surface_model.total_ssite_concentrations = sc; 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.initialize_variables(variables, 0.8, {{"X2", -3}, {"X4", -11}}, std::unordered_map(), 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(); + std::cout << "\t" << sol.free_surface_concentration(0); 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::Vector total_concentrations(raw_data->nb_component()); + total_concentrations << 1, 0.0, 0.6, -2.4, 0.6, 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.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(); + solve_kinetics_problem(); - AutomaticReactionPath test_automatic; + //AutomaticReactionPath test_automatic; - for (int i=0; i<20; ++i) - { - test_automatic.run_step(); - } + //for (int i=0; i<20; ++i) + //{ + // test_automatic.run_step(); + //} return EXIT_SUCCESS; }