diff --git a/tests/specmicp/carboalu.cpp b/tests/specmicp/carboalu.cpp
index 8d574a4..670dce9 100644
--- a/tests/specmicp/carboalu.cpp
+++ b/tests/specmicp/carboalu.cpp
@@ -1,110 +1,115 @@
#include <iostream>
#include "specmicp/reaction_path.hpp"
#include "utils/log.hpp"
double mult = 1.0;
double m_c3s = mult*0.6;
double m_c2s = mult*0.2;
double m_c3a = mult*0.10;
double m_gypsum = mult*0.10;
double wc =0.8;
void solve_carboalu()
{
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
specmicp::database::Database database("data/cemdata_specmicp.js");
std::shared_ptr<specmicp::database::DataContainer> data = database.get_database();
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Al[3+]","Al(OH)4[-]"},
{"Si(OH)4", "SiO(OH)3[-]"}
});
database.swap_components(swapping);
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
double delta_h2co3 = 0.05;
double delta_sio2 = 0.00;
double 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));
std::cout << m_water << std::endl;
model->amount_aqueous = {
{"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water,
delta_h2co3+delta_sio2)},
{"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)},
{"HO[-]", specmicp::reaction_amount_t(0, -delta_h2co3-delta_sio2)},
{"SiO(OH)3[-]", specmicp::reaction_amount_t(0, delta_sio2)},
};
model->amount_minerals = {
{"C3S", specmicp::reaction_amount_t(m_c3s, 0)},
{"C2S", specmicp::reaction_amount_t(m_c2s, 0)},
{"C3A", specmicp::reaction_amount_t(m_c3a, 0)},
{"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0)}
};
model->database_path = "data/cemdata_specmicp.js";
model->minerals_to_keep = {
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2,am",
"Calcite",
"Al(OH)3,am",
"Monosulfoaluminate",
"Tricarboaluminate",
"Monocarboaluminate",
"Hemicarboaluminate",
"Straetlingite",
"Gypsum",
"Ettringite",
"Thaumasite"
};
specmicp::ReactionPathDriver driver(model, data);
driver.dissolve_to_components();
Eigen::VectorXd totaq(data->nb_component);
double totiter = 0;
double totfact = 0;
Eigen::VectorXd x(data->nb_component+data->nb_mineral);
x(0) = m_water;
x.block(1, 0, data->nb_component-1, 1).setConstant(-2);
x.block(data->nb_component, 0, data->nb_mineral, 1).setConstant(0.);
+
//driver.get_options().solver_options.penalization_factor =0.8;
driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
//driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
driver.get_options().solver_options.use_scaling = false;
driver.get_options().solver_options.max_factorization_step = 4;
driver.get_options().solver_options.factor_gradient_search_direction = 100;
driver.get_options().solver_options.max_iter = 50;
driver.get_options().solver_options.maxstep = 50;
driver.get_options().allow_restart = true;
- std::cout << "Reaction_path return_code nb_iter nb_fact pH H2O HO Al C Ca Si SO4";
+ std::cout << "Reaction_path return_code nb_iter nb_fact pH";
+ for (auto it = data->labels_basis.begin(); it!= data->labels_basis.end(); ++it)
+ {
+ std::cout << " " << *it;
+ }
for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it)
{
std::cout << " " << *it;
}
std::cout << std::endl;
const int nb_step = 51;
for (int i=0; i<nb_step; ++i)
{
specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x);
driver.total_aqueous_concentration(x, totaq);
std::cout << i*(delta_sio2+delta_h2co3) << " " << (int) perf.return_code << " "
<< perf.nb_iterations << " " << perf.nb_factorization << " "
<< 14+x(1) << " "
<< x(0) << " " << totaq.block(1,0,data->nb_component-1,1).transpose() << " "
<< x.block(data->nb_component, 0, data->nb_mineral, 1).transpose() << std::endl;
totiter += perf.nb_iterations;
totfact += perf.nb_factorization;
}
std::cout << "Average iterations : " << totiter/nb_step << std::endl;
std::cout << "Average factorization : " << totfact/nb_step << std::endl;
}
int main()
{
for (int i=0; i<1; ++i)
solve_carboalu();
}
diff --git a/tests/specmicp/thermocarbo.cpp b/tests/specmicp/thermocarbo.cpp
index 1fd9696..df69eaa 100644
--- a/tests/specmicp/thermocarbo.cpp
+++ b/tests/specmicp/thermocarbo.cpp
@@ -1,158 +1,171 @@
#include <iostream>
#include "specmicp/reduced_system.hpp"
#include "micpsolver/micpsolver.hpp"
#include "micpsolver/micpsolver_min.hpp"
#include "specmicp/reduced_system_solver.hpp"
#include "database/database.hpp"
#include "specmicp/reaction_path.hpp"
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();
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"}
});
database.swap_components(swapping);
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
double m_c3s = 0.7;
double m_c2s = 0.3;
double wc = 1.0;
double m_water = wc*((3*56.08+60.08)*m_c3s + (2*56.08+60.08)*m_c2s)*1e-3;
double delta_h2co3 = 0.1;
model->amount_aqueous = {
{"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water, delta_h2co3)},
{"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)},
{"HO[-]", specmicp::reaction_amount_t(0, -delta_h2co3)},
};
model->amount_minerals = {
{"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)},
{"C2S", specmicp::reaction_amount_t(m_c2s, 0.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();
+
+
+ Eigen::VectorXd x(data->nb_component+data->nb_mineral);
+ x(0) = 1.0;
+ x.block(1, 0, data->nb_component-1, 1).setConstant(-2.0);
+ x.block(data->nb_component, 0, data->nb_mineral, 1).setZero();
+
//driver.get_options().solver_options.penalization_factor =1.0;
driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
//driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
driver.get_options().solver_options.use_scaling = false;
driver.get_options().solver_options.max_factorization_step = 1;
driver.get_options().solver_options.factor_gradient_search_direction = 50;
driver.get_options().solver_options.maxstep = 50;
driver.get_options().solver_options.maxiter_maxstep = 50;
driver.get_options().solver_options.max_iter = 100;
driver.get_options().allow_restart = true;
- std::cout << "Reaction_path return_code nb_iter nb_fact pH H2O HO C Ca Si "
- << " Portlandite SiO2(am) Jennite Tobermorite Calcite" << std::endl;
- int max_step=41;
+ std::cout << "Reaction_path return_code nb_iter nb_fact pH";
+ for (auto it = data->labels_basis.begin(); it!= data->labels_basis.end(); ++it)
+ {
+ std::cout << " " << *it;
+ }
+ for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it)
+ {
+ std::cout << " " << *it;
+ }
+ std::cout << std::endl;
+ int max_step=41;
for (int i=0; i<max_step; ++i)
{
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 << " " << perf.nb_factorization << " "
<< 14+x(1) << " "
- << x(0) << " " << totaq.block(1,0,4,1).transpose() << " " << x.block(5, 0, 5, 1).transpose() << std::endl;
+ << x(0) << " " << totaq.block(1,0,data->nb_component-1,1).transpose() << " "
+ << x.block(data->nb_component, 0, data->nb_mineral, 1).transpose() << std::endl;
totiter += perf.nb_iterations;
totfact += perf.nb_factorization;
}
std::cout << "Average iterations : " << totiter/max_step << std::endl;
std::cout << "Average factorization : " << totfact/max_step << std::endl;
}
void solve()
{
specmicp::stdlog::ReportLevel() = specmicp::logger::Debug;
specmicp::database::Database database("data/cemdata_specmicp.js");
std::shared_ptr<specmicp::database::DataContainer> data = database.get_database();
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"}
});
database.swap_components(swapping);
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
double m_c3s = 0.7;
double m_c2s = 0.3;
double wc = 1.0;
double m_water = wc*((3*56.08+60.08)*m_c3s + (2*56.08+60.08)*m_c2s)*1e-3;
double delta_h2co3 = 0.1;
model->amount_aqueous = {
{"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water, delta_h2co3)},
{"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)},
{"HO[-]", specmicp::reaction_amount_t(0, -delta_h2co3)},
};
model->amount_minerals = {
{"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)},
{"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)}
};
//x(0) = m_water;
model->database_path = "data/cemdata_specmicp.js";
specmicp::ReactionPathDriver driver(model, data);
driver.dissolve_to_components();
Eigen::VectorXd totaq(data->nb_component);
Eigen::VectorXd x(data->nb_component+data->nb_mineral);
x(0) = 1.0;
x.block(1, 0, data->nb_component-1, 1).setConstant(-2);
- x.block(data->nb_component, 0, data->nb_mineral-1, 1).setZero();
+ x.block(data->nb_component, 0, data->nb_mineral, 1).setZero();
driver.get_options().solver_options.penalization_factor = 0.8;
driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
//driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
driver.get_options().solver_options.use_scaling = false;
driver.get_options().solver_options.max_factorization_step = 1;
driver.get_options().solver_options.factor_gradient_search_direction = 50;
driver.get_options().solver_options.maxstep = 50;
driver.get_options().solver_options.maxiter_maxstep = 50;
driver.get_options().solver_options.max_iter = 100;
driver.get_options().allow_restart = false;
driver.get_options().solver_options.power_descent_condition = 2.0;
- std::cout << "Reaction_path return_code nb_iter nb_fact pH H2O HO C Ca Si "
- << " Portlandite SiO2(am) Jennite Tobermorite Calcite" << std::endl;
+ std::cout << "Reaction_path return_code nb_iter nb_fact pH";
+ for (auto it = data->labels_basis.begin(); it!= data->labels_basis.end(); ++it)
+ {
+ std::cout << " " << *it;
+ }
+ for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it)
+ {
+ std::cout << " " << *it;
+ }
+ std::cout << std::endl;
specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x);
driver.total_aqueous_concentration(x, totaq);
std::cout << 0.0 << " " << (int) perf.return_code << " "
<< perf.nb_iterations << " " << perf.nb_factorization << " "
<< 14+x(1) << " "
<< x(0) << " " << totaq.block(1,0,data->nb_component-1,1).transpose() << " "
<< x.block(data->nb_component, 0, data->nb_mineral, 1).transpose() << std::endl;
std::cout << "Solution \n" << x << std::endl;
}
int main()
{
specmicp::logger::ErrFile::stream() = &std::cerr;
//solve();
solve_lot_reaction_path();
//for (int i=0; i<5000; ++i) solve_lot_reaction_path();
}