Homec4science
Diffusion SpecMiCP / ReactMiCP (kinetics) Browse

specmicp/tests/specmicp/thermocarbo.cpp
d1005cb9e234kinetics

Last Change
Raw File

thermocarbo.cpp
View Options

#include <iostream>
#include "specmicp/specmicp.hpp"
#include "micpsolver/micpsolver.hpp"
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(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::ThermoData> thermoptr = set_thermodata();
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);
//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;
std::cout << "Reaction_path return_code nb_iter pH H2O Ca HO Si C "
<< " 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,
3*m_c3s +2*m_c2s,
5*m_c3s + 3*m_c2s-m_h2co3,
m_c3s + m_c2s,
m_h2co3;
specmicp::micpsolver::MiCPPerformance perf = solve_one(totconc, x, totaq);
std::cout << m_h2co3 << " " << (int) perf.return_code << " " << perf.nb_iterations << " " << 14+x(2) << " "
<< x(0) << " " << totaq.block(1,0,4,1).transpose() << " " << x.block(15, 0, 5, 1).transpose() << std::endl;
}
}
void solve()
{
specmicp::stdlog::ReportLevel() = specmicp::logger::Debug;
std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata();
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,
3*m_c3s +2*m_c2s,
5*m_c3s + 3*m_c2s-m_h2co3,
m_c3s + m_c2s,
m_h2co3;
std::shared_ptr<specmicp::SpecProg> prog = std::make_shared<specmicp::SpecProg>(thermoptr, totconc);
specmicp::micpsolver::MiCPSolver<specmicp::SpecProg> 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;
}
int main()
{
specmicp::logger::ErrFile::stream() = &std::cerr;
//solve();
solve_lot();
//solve_aluminium();
}
c4science · Help