Page MenuHomec4science
Files F67074303

alkaliactivated.cpp
No OneTemporary
Actions

File Metadata

Created
Wed, Jun 19, 18:36

alkaliactivated.cpp
View Options

#include <iostream>
#include "specmicp/reaction_path.hpp"
#include "utils/log.hpp"
void solve_lot_reaction_path()
{
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 mult = 2;
double m_naoh = 0.1;
double m_sio2 = mult*0.5;
double m_al2o3 = mult*0.5;
double m_cao = mult*1.0;
double m_gypsum = 0.2;
//double wc = 0.5;
double m_water = 1.0;
//double delta_h2co3 = 0.1;
double delta_cao = 0.1;
model->amount_aqueous = {
{"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water -2*m_al2o3-m_cao-m_sio2, +delta_cao)},
//{"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)},
{"HO[-]", specmicp::reaction_amount_t(m_naoh - 2*m_al2o3 +2*m_cao-m_sio2, -2*delta_cao)},
{"Na[+]", specmicp::reaction_amount_t(m_naoh, 0)},
{"Al(OH)4[-]", specmicp::reaction_amount_t(2*m_al2o3, 0)},
{"SiO(OH)3[-]", specmicp::reaction_amount_t(m_sio2, 0)},
{"Ca[2+]", specmicp::reaction_amount_t(m_cao, -delta_cao)}
};
model->amount_minerals = {
{"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0)}
};
model->database_path = "data/cemdata_specmicp.js";
double totiter = 0;
double totfact = 0;
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();
//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 = 2;
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";
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=11;
for (int i=0; i<max_step; ++i)
{
specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x);
specmicp::EquilibriumState solution = driver.get_current_solution();
solution.total_aqueous_concentrations(totaq);
std::cout << i*delta_cao << " " << (int) perf.return_code << " "
<< perf.nb_iterations << " " << perf.nb_factorization << " "
<< solution.pH() << " "
<< solution.mass_water() << " " << totaq.block(1,0,data->nb_component-1,1).transpose();
for (int m=0; m < data->nb_mineral; ++m) std::cout << " " << solution.moles_mineral(m);
std::cout << 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;
}
int main()
{
solve_lot_reaction_path();
}

Event Timeline

c4science · Help