Page MenuHomec4science
Files F68703836

thermocarbo.cpp
No OneTemporary
Actions

File Metadata

Created
Fri, Jun 28, 14:14

thermocarbo.cpp
View Options

#include <iostream>
#include "specmicp/extended_system.hpp"
#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) = -6.0;
x(2) = -6.0;
x(3) = -6.0;
x(4) = -6.0;
x(5) = 1.0;
x(6) = 1.0;
x(7) = 1.0;
x(8) = 1.0;
x(9) = 1.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();
//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;
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;
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;
Eigen::VectorXd x(10);
x(0) = 1;
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;
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";
Eigen::VectorXd totaq(5);
specmicp::ReactionPathDriver driver(model, data);
driver.dissolve_to_components();
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;
specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x);
driver.total_aqueous_concentration(x, totaq);
std::cout << 0*(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;
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();
}

Event Timeline

c4science · Help