diff --git a/src/specmicp/reaction_path.cpp b/src/specmicp/reaction_path.cpp
index 7f669cb..00c8a03 100644
--- a/src/specmicp/reaction_path.cpp
+++ b/src/specmicp/reaction_path.cpp
@@ -1,83 +1,82 @@
/*-------------------------------------------------------
- Module : specmicp
- File : reaction_path.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "reaction_path.hpp"
#include "thermodata.hpp"
#include "micpsolver/micpsolver.hpp"
#include "utils/log.hpp"
#include <iostream>
namespace specmicp {
void ReactionPathDriver::read_database()
{
m_database.parse_database(m_model->database_path);
m_database.make_canonical();
m_data = m_database.get_database();
}
int ReactionPathDriver::label_to_id(const std::string& label, const std::vector<std::string>& list_labels)
{
auto it = std::find(list_labels.begin(), list_labels.end(), label);
if (it == list_labels.end())
{
throw std::invalid_argument("Unknow species : " + label);
}
return it - list_labels.begin();
}
void ReactionPathDriver::dissolve_to_components()
{
m_tot_conc = Eigen::VectorXd(m_data->nb_component);
m_tot_conc_increase = Eigen::VectorXd(m_data->nb_component);
// everything is dissolved ...
for (auto it=m_model->amount_components.begin(); it!=m_model->amount_components.end(); ++it)
{
const int id = label_to_id(it->first, m_data->labels_basis);
m_tot_conc(id) = it->second.first;
m_tot_conc_increase(id) = it->second.second;
}
std::vector<int> to_remove;
to_remove.reserve(10);
int new_i = 0;
for (int i=0; i<m_data->nb_component; ++i)
{
if (m_tot_conc(i) == 0 and m_tot_conc_increase(i) == 0)
{
to_remove.push_back(0);
continue;
}
m_tot_conc(new_i) = m_tot_conc(i);
m_tot_conc_increase(new_i) = m_tot_conc_increase(i);
++new_i;
}
m_tot_conc.conservativeResize(new_i);
m_tot_conc_increase.conservativeResize(new_i);
- std::cout << m_tot_conc << std::endl;
- std::cout << m_tot_conc_increase << std::endl;
m_database.remove_components(to_remove);
}
micpsolver::MiCPPerformance ReactionPathDriver::one_step(Eigen::VectorXd& x)
{
- m_current_solver = ReducedSystemSolver (m_data, m_tot_conc, m_options);
+ m_current_solver = ReducedSystemSolver(m_data, m_tot_conc);
+ m_current_solver.get_options().solver_options = m_solver_options;
micpsolver::MiCPPerformance perf = m_current_solver.solve(x);
if (perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized)
{
ERROR << "Failed to solve the system !";
}
m_tot_conc += m_tot_conc_increase;
return perf;
}
} // end namespace specmicp
diff --git a/src/specmicp/reaction_path.hpp b/src/specmicp/reaction_path.hpp
index a5cdfec..c78df8c 100644
--- a/src/specmicp/reaction_path.hpp
+++ b/src/specmicp/reaction_path.hpp
@@ -1,106 +1,106 @@
/*-------------------------------------------------------
- Module : specmicp
- File : reaction_path.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMICP_SPECMICP_REACTIONPATH_HPP
#define SPECMICP_SPECMICP_REACTIONPATH_HPP
//! \file reaction_path.hpp reaction path driver
#include <memory>
#include <map>
#include "database/database.hpp"
#include "micpsolver/micpsolver_structs.hpp"
#include "reduced_system_solver.hpp"
namespace specmicp {
//! \brief typedef to a reaction path amount
//!
//! The first element is the inital amount,
//! The second element is the increase to add at every step
using reaction_amount_t = std::pair<double, double>;
//! \brief A reaction path model
//!
//! Contains every data we need to run a reaction path model
struct ReactionPathModel
{
//! The amount of components
std::map<std::string, reaction_amount_t> amount_components;
//! The amount of aqueous species
std::map<std::string, reaction_amount_t> amount_aqueous;
//! The amount of minerals
std::map<std::string, reaction_amount_t> amount_minerals;
//! the number of step
int nb_step;
//! the path to the database - optional
std::string database_path;
};
//! \brief Driver to solve a reaction path model
class ReactionPathDriver
{
public:
//! \brief Initialise from a model
//!
//! Build a database from the attribute database_path of the model
ReactionPathDriver(std::shared_ptr<ReactionPathModel> model):
m_model(model)
{
read_database();
}
//! \brief Initialise the driver from a model and a raw database
ReactionPathDriver(std::shared_ptr<ReactionPathModel> model,
std::shared_ptr<database::DataContainer> data):
m_model(model),
m_database(data),
m_data(data)
{}
//! Initialize the database
void read_database();
//! Dissolve everything into components
void dissolve_to_components();
//! \brief Return the id of a species in 'list_labels' from its id
int label_to_id(const std::string& label, const std::vector<std::string>& list_labels);
//! \brief Perform one step
micpsolver::MiCPPerformance one_step(Eigen::VectorXd& x);
//! Set total aqueous concentrations
void total_aqueous_concentration(const Eigen::VectorXd& x, Eigen::VectorXd& totaq) {
m_current_solver.total_aqueous_concentration(x, totaq);
}
private:
std::shared_ptr<ReactionPathModel> m_model;
specmicp::database::Database m_database;
std::shared_ptr<specmicp::database::DataContainer> m_data;
Eigen::VectorXd m_tot_conc;
Eigen::VectorXd m_tot_conc_increase;
- micpsolver::MiCPSolverOptions m_options;
+ micpsolver::MiCPSolverOptions m_solver_options;
ReducedSystemSolver m_current_solver;
//std::shared_ptr<ThermoData> m_thermo;
};
} // end namespace specmicp
#endif // SPECMICP_SPECMICP_REACTIONPATH_HPP
diff --git a/src/specmicp/reduced_system_solver.cpp b/src/specmicp/reduced_system_solver.cpp
index dc2d03d..70e6828 100644
--- a/src/specmicp/reduced_system_solver.cpp
+++ b/src/specmicp/reduced_system_solver.cpp
@@ -1,122 +1,130 @@
/*-------------------------------------------------------
- Module : specmicp
- File : reduced_system_solver
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "reduced_system_solver.hpp"
#include "micpsolver/micpsolver.hpp"
#include <iostream>
namespace specmicp {
-ReducedSystemSolver::ReducedSystemSolver(
- std::shared_ptr<database::DataContainer> data,
- Eigen::VectorXd& totconc,
- micpsolver::MiCPSolverOptions options):
+ReducedSystemSolver::ReducedSystemSolver(std::shared_ptr<database::DataContainer> data,
+ Eigen::VectorXd& totconc):
m_data(data),
m_system(std::make_shared<ReducedSystem>(std::make_shared<ThermoData>(data), totconc)),
m_var(Eigen::VectorXd::Zero(data->nb_component+data->nb_mineral)),
- m_options(options)
+ m_options()
{}
micpsolver::MiCPPerformance ReducedSystemSolver::solve(Eigen::VectorXd &x)
{
set_true_variable_vector(x);
micpsolver::MiCPPerformance perf = solve();
set_return_vector(x);
return perf;
}
void ReducedSystemSolver::set_true_variable_vector(const Eigen::VectorXd &x)
{
const std::vector<int>& non_active_component = m_system->get_non_active_component();
if (non_active_component.size() == 0)
{
m_var = x;
for (int i=0; i<m_var.rows(); ++i)
{
if (m_var(i) == -HUGE_VAL)
{
m_var(i) = -6;
}
}
}
else
{
unsigned int new_i = 0;
for (int i=0; i<m_data->nb_component; ++i)
{
auto it = std::find(non_active_component.begin(), non_active_component.end(),i);
if (it != non_active_component.end()) continue;
m_var(new_i) = x(i);
++new_i;
}
assert(new_i == m_data->nb_component - non_active_component.size());
for (int m=0; m<m_data->nb_mineral; ++m)
{
for (auto it=non_active_component.begin(); it!=non_active_component.end(); ++it)
{
if (m_data->nu_mineral(m, *it) != 0) continue;
m_var(new_i) = x(m_data->nb_component+m);
++new_i;
}
}
m_var.conservativeResize(new_i);
}
}
micpsolver::MiCPPerformance ReducedSystemSolver::solve()
{
- micpsolver::MiCPSolver<ReducedSystem> solver(m_system);
- solver.set_options(m_options);
- solver.solve(m_var);
- return solver.get_performance();
+ if (get_options().ncp_function == micpsolver::NCPfunction::penalizedFB)
+ {
+ micpsolver::MiCPSolver<ReducedSystem,micpsolver::NCPfunction::penalizedFB> solver(m_system);
+ solver.set_options(get_options().solver_options);
+ solver.solve(m_var);
+ return solver.get_performance();
+ }
+ else
+ {
+ micpsolver::MiCPSolver<ReducedSystem,micpsolver::NCPfunction::min> solver(m_system);
+ solver.set_options(get_options().solver_options);
+ solver.solve(m_var);
+ return solver.get_performance();
+ }
}
void ReducedSystemSolver::set_return_vector(Eigen::VectorXd &x)
{
const std::vector<int>& non_active_component = m_system->get_non_active_component();
if (non_active_component.size() == 0)
{
x = m_var;
}
else
{
unsigned int new_i = 0;
for (int i=0; i<m_data->nb_component; ++i)
{
auto it = std::find(non_active_component.begin(), non_active_component.end(),i);
if (it != non_active_component.end())
{
x(i) = -HUGE_VAL; // FIXME
continue;
}
x(i) = m_var(new_i) ;
++new_i;
}
assert(new_i == m_data->nb_component - non_active_component.size());
for (int m=0; m<m_data->nb_mineral; ++m)
{
for (auto it=non_active_component.begin(); it!=non_active_component.end(); ++it)
{
if (m_data->nu_mineral(m, *it) != 0)
{
x(m_data->nb_component+m) = 0;
continue;
}
x(m_data->nb_component+m) =m_var(new_i);
++new_i;
}
}
}
}
} // end namespace specmicp
diff --git a/src/specmicp/reduced_system_solver.hpp b/src/specmicp/reduced_system_solver.hpp
index 7f49f25..593af28 100644
--- a/src/specmicp/reduced_system_solver.hpp
+++ b/src/specmicp/reduced_system_solver.hpp
@@ -1,76 +1,88 @@
/*-------------------------------------------------------
- Module : specmicp
- File : reduced_system_solver.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMICP_SPECMICP_REDUCEDSYSTEMSOLVER_HPP
#define SPECMICP_SPECMICP_REDUCEDSYSTEMSOLVER_HPP
#include "database/data_container.hpp"
#include "reduced_system.hpp"
#include "micpsolver/micpsolver_structs.hpp"
//! \file reduced_system_solver.hpp Solve a reduced system
namespace specmicp {
+struct ReducedSystemSolverOptions
+{
+ micpsolver::NCPfunction ncp_function;
+ micpsolver::MiCPSolverOptions solver_options;
+
+ ReducedSystemSolverOptions():
+ ncp_function(micpsolver::NCPfunction::penalizedFB),
+ solver_options()
+ {}
+};
//! \brief Solver a reduced system
//!
//! Take care of possibly non-existing component in the system
class ReducedSystemSolver
{
public:
//! Default constructor - you need to use another method to initialize it correctly
ReducedSystemSolver() {}
ReducedSystemSolver(std::shared_ptr<database::DataContainer> data,
- Eigen::VectorXd& totconc,
- micpsolver::MiCPSolverOptions options=micpsolver::MiCPSolverOptions());
+ Eigen::VectorXd& totconc);
//! \brief solve the problem using initial guess x
//!
//! @param[in,out] x, in -> initial guess, out -> solution
micpsolver::MiCPPerformance solve(Eigen::VectorXd& x);
//! Return the system used for the computation
std::shared_ptr<ReducedSystem> get_system() {return m_system;}
//! Set total aqueous concentrations
void total_aqueous_concentration(const Eigen::VectorXd& x, Eigen::VectorXd& totaq)
{
m_system->aqueous_tot_conc(x, totaq);
}
+ const ReducedSystemSolverOptions& get_options() const {return m_options;}
+ ReducedSystemSolverOptions& get_options() {return m_options;}
+
private:
//! \brief set up the true variable vector
void set_true_variable_vector(const Eigen::VectorXd& x);
//! \brief set up the true solution vector
//!
//! add zero components
void set_return_vector(Eigen::VectorXd& x);
//! \brief solve the problem
micpsolver::MiCPPerformance solve();
+ ReducedSystemSolverOptions m_options;
std::shared_ptr<database::DataContainer> m_data;
std::shared_ptr<ReducedSystem> m_system;
Eigen::VectorXd m_var;
- micpsolver::MiCPSolverOptions m_options;
};
inline micpsolver::MiCPPerformance solve_reduced_system(std::shared_ptr<database::DataContainer> data,
Eigen::VectorXd& totconc,
Eigen::VectorXd& x)
{
return ReducedSystemSolver(data, totconc).solve(x);
}
} // end namespace specmicp
#endif // SPECMICP_SPECMICP_REDUCEDSYSTEMSOLVER_HPP
diff --git a/tests/specmicp/thermocarbo.cpp b/tests/specmicp/thermocarbo.cpp
index 5a95507..2d7f446 100644
--- a/tests/specmicp/thermocarbo.cpp
+++ b/tests/specmicp/thermocarbo.cpp
@@ -1,469 +1,533 @@
#include <iostream>
#include "specmicp/extended_system.hpp"
#include "specmicp/reduced_system.hpp"
#include "micpsolver/micpsolver.hpp"
#include "specmicp/reduced_system_solver.hpp"
#include "database/database.hpp"
#include "specmicp/reaction_path.hpp"
std::shared_ptr<specmicp::ThermoData> set_thermodata_from_database()
{
specmicp::database::Database database("data/cemdata_specmicp.js");
std::shared_ptr<specmicp::database::DataContainer> data = database.get_database();
//std::cout << data->nu_aqueous << std::endl;
//std::cout << data->nu_mineral << std::endl;
database.change_basis(std::vector<int>({0, 5, 2, 3, 8}));
//std::cout << data->nu_aqueous << std::endl;
//std::cout << data->nu_mineral << std::endl;
std::cout << "Basis : " << std::endl;
for (auto it = data->labels_basis.begin(); it!=data->labels_basis.end(); ++it)
{
std::cout << *it << "\t" << data->param_aq.row(it - data->labels_basis.begin()) << std::endl;
}
std::cout << std::endl;
std::cout << "Aqueous : \n";
for (auto it = data->labels_aqueous.begin(); it!=data->labels_aqueous.end(); ++it)
{
const int j = it - data->labels_aqueous.begin();
std::cout << *it << "\t\t" << data->nu_aqueous.row(j) << "\t"
<< data->logk_aqueous(j) << "\t" << data->param_aq.row(j+data->nb_component) << std::endl;
}
std::cout << "Minerals : \n";
for (auto it = data->labels_minerals.begin(); it!=data->labels_minerals.end(); ++it)
{
const int j = it - data->labels_minerals.begin();
std::cout << *it << "\t\t" << data->nu_mineral.row(j) << "\t" << data->logk_mineral(j) << std::endl;
}
std::cout << std::endl;
std::shared_ptr<specmicp::ThermoData> ptr = std::make_shared<specmicp::ThermoData>(data);
return ptr;
}
//std::shared_ptr<specmicp::ThermoData> set_thermodata()
//{
// int nmin = 5;
// Eigen::MatrixXd nu(10+nmin, 5);
// nu << 1, 0, -1, 0, 0,
// 1, 0, -1, 1, 0,
// -1, 0, 1, 1, 0,
// 0, 1, 1, 0, 0,
// -1, 0, 1, 0, 1,
// 0, 0, -1, 0, 1,
// -1, 1, 1, 0, 1,
// 0, 1, 0, 0, 1,
// 0, 1, 0, 1, 0,
// -1, 1, 1, 1, 0,
// // minerals
// 0, 1, 2, 0, 0,
// -1, 0, -1, 1, 0,
// -0.5567, 1.6667, 2.3333, 1, 0,
// -0.5, 0.8333, 0.6667, 1, 0,
// -1, 1, 1, 0, 1
// ;
// Eigen::VectorXd logkr(10+nmin);
// logkr << 13.9995, 4.8595, 0.1005, -1.2195,
// -3.6706, 7.6476, -6.8947, -1.1057,
// -1.2, -4.4995,
// // minerals
// -5.1995,
// 1.476,
// -13.17,
// -8.0,
// -12.1505;
// Eigen::MatrixXd param_aqueous(Eigen::MatrixXd::Zero(5+10,3));
// param_aqueous << 0, 0.0 , 0.0 , // H2O
// 2, 4.86, 0.15, // Ca+2
// -1, 10.65, 0.0 , // HO-
// -1, 4.0 , 0.0 , // SiO(HO)3-
// -1, 5.40, 0.0 , // HCO3-
// 1, 9.0 , 0.0 , // H+
// 0, 0.0 , 0.1 , // Si(OH)4
// -2, 4.0 , 0.0 , // SiO2(OH)-2
// 1, 4.0 , 0.0 , // CaOH+
// -2, 5.40, 0.0 , // CO3-2
// 0, 0.0 , 0.1 , // CO2(aq)
// 0, 0.0 , 0.1 , // CaCO3(aq)
// 1, 4.0 , 0.0 , // CaHCO3+
// 1, 4.0 , 0.0 , // CaSiO(HO)3+
// 0, 0.0 , 0.1 // CaSiO2(HO)2(aq)
// ;
// std::shared_ptr<specmicp::ThermoData> ptr = std::make_shared<specmicp::ThermoData>(10, nmin, nu, logkr, param_aqueous);
// return ptr;
//}
//std::shared_ptr<specmicp::ThermoData> set_thermodata_aluminum()
//{
// Eigen::MatrixXd nu(5, 3);
// nu << 1, 1, -1,
// 2, 1, -2,
// 3, 1, -3,
// 4, 1, -4,
// 1, 0, -1
// ;
// Eigen::VectorXd logkr(5);
// logkr << 5.0, 10.1, 16.9, 22.7, 14.0;
// Eigen::MatrixXd paramaq(Eigen::MatrixXd::Zero(5+3,3)); // ideal case
// std::shared_ptr<specmicp::ThermoData> ptr = std::make_shared<specmicp::ThermoData>(5, 0, nu, logkr, paramaq);
// return ptr;
//}
//void solve_aluminium()
//{
// std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata_aluminum();
// Eigen::VectorXd totconc(3);
// totconc << 1.0/specmicp::molar_mass_water, 1e-6, -3e-6;
// //totconc << 1e-6, -3e-6;
// std::shared_ptr<specmicp::SpecProg> prog = std::make_shared<specmicp::SpecProg>(thermoptr, totconc, true);
// specmicp::micpsolver::MiCPSolver<specmicp::SpecProg> solver(prog);
// //solver.get_options().fvectol = 1e-6;
// solver.get_options().maxstep = 20.0;
// solver.get_options().maxiter_maxstep = 100;
// //solver.get_options().factor_descent_condition = -1e-6;
// Eigen::VectorXd x(prog->total_variables());
// //x.setConstant(-1);
// x(0) = 1.0;
// x(1) = -9;
// x(2) = -13;
// prog->init_secondary_conc(x);
// std::cout << x << std::endl;
//// Eigen::VectorXd residual(20);
//// prog->get_residuals(x, residual);
//// std::cout << "Residual \n ------ \n" << residual << std::endl;
//// Eigen::MatrixXd jacob(20, 20);
//// prog->get_jacobian(x, jacob);
//// std::cout << "Jacobian \n ------ \n" << jacob << std::endl;
// std::cout << "Return Code : " << (int) solver.solve(x) << std::endl;
// std::cout << "Solution \n -------- \n" << x << std::endl;
//}
specmicp::micpsolver::MiCPPerformance solve_one(
std::shared_ptr<specmicp::ThermoData> thermoptr,
const Eigen::VectorXd& totconc,
Eigen::VectorXd& x,
Eigen::VectorXd& totaq)
{
x(15) = std::max(0.0, x(15)-1.0);
x(16) = std::max(0.0, x(16)-1.0);
x(17) = std::max(0.0, x(17)-1.0);
x(18) = std::max(0.0, x(18)-1.0);
x(19) = std::max(0.0, x(19)-1.0);
std::shared_ptr<specmicp::SpecProg> prog = std::make_shared<specmicp::SpecProg>(thermoptr, totconc);
prog->init_secondary_conc(x);
- specmicp::micpsolver::MiCPSolver<specmicp::SpecProg> solver(prog);
+ specmicp::micpsolver::MiCPSolver<specmicp::SpecProg, specmicp::micpsolver::NCPfunction::penalizedFB> solver(prog);
//solver.get_options().fvectol = 1e-6;
solver.get_options().maxstep = 200.0;
solver.get_options().maxiter_maxstep = 100;
solver.get_options().use_crashing = false;
solver.get_options().use_scaling = false;
solver.get_options().projection_min_variable = 1e-6;
solver.solve(x);
prog->aqueous_tot_conc(x, totaq);
return solver.get_performance();
}
void solve_lot()
{
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
Eigen::VectorXd x(20);
x(0) = 10.0;
x(1) = -2.0;
x(2) = -2.0;
x(3) = -2.0;
x(4) = -2.0;
x(15) = 0.0;
x(16) = 0.0;
x(17) = 0.0;
x(18) = 0.0;
x(19) = 0.0;
//x.setConstant(-1);
//prog->init_secondary_conc(x);
Eigen::VectorXd totaq(5);
double m_c3s = 0.70;
double m_c2s = 0.3;
double m_h2co3 = 0.0;
double totiter = 0;
double totfact = 0;
std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata_from_database();
std::cout << "Reaction_path return_code nb_iter pH H2O HO C Ca Si "
<< " Portlandite SiO2(am) Jennite Tobermorite Calcite" << std::endl;
for (int i=0; i<40; ++i)
{
m_h2co3 = 0.1*(1+i);
Eigen::VectorXd totconc(5);
totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
5*m_c3s + 3*m_c2s-m_h2co3,
m_h2co3,
3*m_c3s +2*m_c2s,
m_c3s + m_c2s;
specmicp::micpsolver::MiCPPerformance perf = solve_one(thermoptr, totconc, x, totaq);
std::cout << m_h2co3 << " " << (int) perf.return_code << " " << perf.nb_iterations << " " << 14+x(1) << " "
<< x(0) << " " << totaq.block(1,0,4,1).transpose() << " " << x.block(15, 0, 5, 1).transpose() << std::endl;
totiter += perf.nb_iterations;
totfact += perf.nb_factorization;
}
std::cout << "Average iterations : " << totiter/40 << std::endl;
std::cout << "Average factorization : " << totfact/40 << std::endl;
}
specmicp::micpsolver::MiCPPerformance solve_one_reduced(
std::shared_ptr<specmicp::database::DataContainer> data,
Eigen::VectorXd& totconc,
Eigen::VectorXd& x,
Eigen::VectorXd& totaq)
{
x(5) = std::max(0.0, x(5)-1.0);
x(6) = std::max(0.0, x(6)-1.0);
x(7) = std::max(0.0, x(7)-1.0);
x(8) = std::max(0.0, x(8)-1.0);
x(9) = std::max(0.0, x(9)-1.0);
specmicp::micpsolver::MiCPSolverOptions options;
options.maxstep = 100.0;
options.maxiter_maxstep = 100;
options.use_crashing = false;
options.use_scaling = true;
options.projection_min_variable = 1e-6;
- specmicp::ReducedSystemSolver solver(data, totconc, options);
+ specmicp::ReducedSystemSolver solver(data, totconc);
+ solver.get_options().solver_options = options;
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x);
solver.total_aqueous_concentration(x, totaq);
return perf;
}
void solve_lot_reduced()
{
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
Eigen::VectorXd x(10);
x(0) = 1.0;
x(1) = -2.0;
x(2) = -2.0;
x(3) = -2.0;
x(4) = -2.0;
x(5) = 0.0;
x(6) = 0.0;
x(7) = 0.0;
x(8) = 0.0;
x(9) = 0.0;
//x.setConstant(-1);
//prog->init_secondary_conc(x);
Eigen::VectorXd totaq(5);
double m_c3s = 0.70;
double m_c2s = 0.3;
double m_h2co3 = 0.0;
double totiter = 0;
double totfact = 0;
std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata_from_database();
Eigen::VectorXd totconc(5);
totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
5*m_c3s + 3*m_c2s-m_h2co3,
m_h2co3,
3*m_c3s +2*m_c2s,
m_c3s + m_c2s;
std::cout << "Reaction_path return_code nb_iter pH H2O HO C Ca Si "
<< " Portlandite SiO2(am) Jennite Tobermorite Calcite" << std::endl;
for (int i=0; i<41; ++i)
{
m_h2co3 = 0.1*(i);
totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
5*m_c3s + 3*m_c2s-m_h2co3,
m_h2co3,
3*m_c3s +2*m_c2s,
m_c3s + m_c2s;
specmicp::micpsolver::MiCPPerformance perf = solve_one_reduced(thermoptr->get_raw_database(), totconc, x, totaq);
//std::cout << x << std::endl;
std::cout << m_h2co3 << " " << (int) perf.return_code << " " << perf.nb_iterations << " " << 14+x(1) << " "
<< x(0) << " " << totaq.block(1,0,4,1).transpose() << " " << x.block(5, 0, 5, 1).transpose() << std::endl;
totiter += perf.nb_iterations;
totfact += perf.nb_factorization;
}
std::cout << "Average iterations : " << totiter/41 << std::endl;
std::cout << "Average factorization : " << totfact/41 << std::endl;
}
void solve_lot_reaction_path()
{ specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
Eigen::VectorXd x(10);
x(0) = 1.0;
x(1) = -2.0;
x(2) = -2.0;
x(3) = -2.0;
x(4) = -2.0;
x(5) = 0.0;
x(6) = 0.0;
x(7) = 0.0;
x(8) = 0.0;
x(9) = 0.0;
specmicp::database::Database database("data/cemdata_specmicp.js");
std::shared_ptr<specmicp::database::DataContainer> data = database.get_database();
database.change_basis(std::vector<int>({0, 5, 2, 3, 8}));
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
double m_c3s = 0.70;
double m_c2s = 0.3;
double delta_h2co3 = 0.1;
model->amount_components = {
{"H2O", specmicp::reaction_amount_t(1/specmicp::molar_mass_water- 4*m_c3s - 3*m_c2s, delta_h2co3)},
{"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)},
{"Ca[2+]", specmicp::reaction_amount_t(3*m_c3s +2*m_c2s, 0)},
{"HO[-]", specmicp::reaction_amount_t(5*m_c3s + 3*m_c2s, -delta_h2co3)},
{"SiO(OH)3[-]", specmicp::reaction_amount_t(m_c3s + m_c2s, 0)}
};
model->database_path = "data/cemdata_specmicp.js";
Eigen::VectorXd totaq(5);
double totiter = 0;
double totfact = 0;
specmicp::ReactionPathDriver driver(model, data);
driver.dissolve_to_components();
std::cout << "Reaction_path return_code nb_iter pH H2O HO C Ca Si "
<< " Portlandite SiO2(am) Jennite Tobermorite Calcite" << std::endl;
for (int i=0; i<41; ++i)
{
+ x(0) = 1.0;
+ x(1) = -2.0;
+ x(2) = -2.0;
+ x(3) = -2.0;
+ x(4) = -2.0;
+ x(5) = 0.0;
+ x(6) = 0.0;
+ x(7) = 0.0;
+ x(8) = 0.0;
+ x(9) = 0.0;
specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x);
driver.total_aqueous_concentration(x, totaq);
std::cout << i*(0.1) << " " << (int) perf.return_code << " " << perf.nb_iterations << " " << 14+x(1) << " "
<< x(0) << " " << totaq.block(1,0,4,1).transpose() << " " << x.block(5, 0, 5, 1).transpose() << std::endl;
totiter += perf.nb_iterations;
totfact += perf.nb_factorization;
}
std::cout << "Average iterations : " << totiter/41 << std::endl;
std::cout << "Average factorization : " << totfact/41 << std::endl;
}
void solve()
{
specmicp::stdlog::ReportLevel() = specmicp::logger::Debug;
std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata_from_database();
Eigen::VectorXd totconc(5);
double m_c3s = 0.70;
double m_c2s = 0.3;
double m_h2co3 = 1.9;
totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
5*m_c3s + 3*m_c2s-m_h2co3,
m_h2co3,
3*m_c3s +2*m_c2s,
m_c3s + m_c2s;
std::shared_ptr<specmicp::SpecProg> prog = std::make_shared<specmicp::SpecProg>(thermoptr, totconc);
- specmicp::micpsolver::MiCPSolver<specmicp::SpecProg> solver(prog);
+ specmicp::micpsolver::MiCPSolver<specmicp::SpecProg, specmicp::micpsolver::NCPfunction::penalizedFB> solver(prog);
//solver.get_options().fvectol = 1e-6;
solver.get_options().max_iter = 100;
solver.get_options().maxstep = 200.0;
solver.get_options().maxiter_maxstep = 100;
solver.get_options().factor_descent_condition = 1e-4;
solver.get_options().use_crashing = 0;
Eigen::VectorXd x(prog->total_variables());
//x.setConstant(-1);
x(0) = 1.0;
x(1) = -2.0;
x(2) = -2.0;
x(3) = -2.0;
x(4) = -2.0;
x(15) = 0.0;
x(16) = 0.0;
x(17) = 0.0;
x(18) = 0.0;
x(19) = 0.0;
prog->init_secondary_conc(x);
std::cout << x << std::endl;
// Eigen::VectorXd residual(20);
// prog->get_residuals(x, residual);
// std::cout << "Residual \n ------ \n" << residual << std::endl;
// Eigen::MatrixXd jacob(20, 20);
// prog->get_jacobian(x, jacob);
// std::cout << "Jacobian \n ------ \n" << jacob << std::endl;
std::cout << "Return Code : " << (int) solver.solve(x) << std::endl;
std::cout << "Solution \n -------- \n" << x << std::endl;
}
void solve_reduced()
{
specmicp::stdlog::ReportLevel() = specmicp::logger::Debug;
std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata_from_database();
Eigen::VectorXd totconc(5);
double m_c3s = 0.70;
double m_c2s = 0.3;
double m_h2co3 = 2.8;
totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
5*m_c3s + 3*m_c2s-m_h2co3,
m_h2co3,
3*m_c3s +2*m_c2s,
m_c3s + m_c2s;
std::shared_ptr<specmicp::ReducedSystem> prog = std::make_shared<specmicp::ReducedSystem>(thermoptr, totconc);
- specmicp::micpsolver::MiCPSolver<specmicp::ReducedSystem> solver(prog);
+ specmicp::micpsolver::MiCPSolver<specmicp::ReducedSystem, specmicp::micpsolver::NCPfunction::penalizedFB> solver(prog);
//solver.get_options().fvectol = 1e-6;
solver.get_options().max_iter = 100;
solver.get_options().maxstep = 200;
solver.get_options().maxiter_maxstep = 100;
solver.get_options().factor_descent_condition = 1e-4;
solver.get_options().use_crashing = false;
solver.get_options().use_scaling = true;
Eigen::VectorXd x(prog->total_variables());
//x.setConstant(-1);
x(0) = 1.0485;
x(1) = -5.01517;
x(2) = -3.52536;
x(3) = -3.50741;
x(4) = -3.5383;
x(5) = 0.0;
x(6) = 0.997242;
x(7) = 0.0;
x(8) = 0.0;
x(9) = 2.69967;
std::cout << x << std::endl;
// Eigen::VectorXd residual(20);
// prog->get_residuals(x, residual);
// std::cout << "Residual \n ------ \n" << residual << std::endl;
// Eigen::MatrixXd jacob(20, 20);
// prog->get_jacobian(x, jacob);
// std::cout << "Jacobian \n ------ \n" << jacob << std::endl;
std::cout << "Return Code : " << (int) solver.solve(x) << std::endl;
std::cout << "Solution \n -------- \n" << x << std::endl;
}
+
+void solve_reduced_min()
+{
+
+ specmicp::stdlog::ReportLevel() = specmicp::logger::Debug;
+ std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata_from_database();
+ Eigen::VectorXd totconc(5);
+ double m_c3s = 0.70;
+ double m_c2s = 0.3;
+ double m_h2co3 = 2.8;
+ totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
+ 5*m_c3s + 3*m_c2s-m_h2co3,
+ m_h2co3,
+ 3*m_c3s +2*m_c2s,
+ m_c3s + m_c2s;
+ std::shared_ptr<specmicp::ReducedSystem> prog = std::make_shared<specmicp::ReducedSystem>(thermoptr, totconc);
+ specmicp::micpsolver::MiCPSolver<specmicp::ReducedSystem, specmicp::micpsolver::NCPfunction::min> solver(prog);
+ //solver.get_options().fvectol = 1e-6;
+ solver.get_options().max_iter = 100;
+ solver.get_options().maxstep = 200;
+ solver.get_options().maxiter_maxstep = 100;
+ solver.get_options().factor_descent_condition = 1e-4;
+ solver.get_options().use_crashing = false;
+ solver.get_options().use_scaling = true;
+ Eigen::VectorXd x(prog->total_variables());
+ //x.setConstant(-1);
+ x(0) = 1.0;
+ x(1) = -2;
+ x(2) = -2;
+ x(3) = -2;
+ x(4) = -2;
+ x(5) = 0.0;
+ x(6) = 0.0;
+ x(7) = 0.0;
+ x(8) = 0.0;
+ x(9) = 0.0;
+
+ std::cout << x << std::endl;
+
+// Eigen::VectorXd residual(20);
+// prog->get_residuals(x, residual);
+// std::cout << "Residual \n ------ \n" << residual << std::endl;
+// Eigen::MatrixXd jacob(20, 20);
+// prog->get_jacobian(x, jacob);
+// std::cout << "Jacobian \n ------ \n" << jacob << std::endl;
+
+ std::cout << "Return Code : " << (int) solver.solve(x) << std::endl;
+
+
+ std::cout << "Solution \n -------- \n" << x << std::endl;
+}
+
int main()
{
specmicp::logger::ErrFile::stream() = &std::cerr;
//set_thermodata_from_database();
//solve();
//solve_lot();
//solve_reduced();
+ //solve_reduced_min();
//solve_lot_reduced();
//solve_aluminium();
solve_lot_reaction_path();
}