Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/tests/specmicp/thermocarbo.cpp
1c6661f02cc4master

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.0 , // 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.0 , // CO2(aq)
0, 0.0 , 0.0 , // CaCO3(aq)
1, 0.0 , 0.0 , // CaHCO3+
1, 4.0 , 0.0 , // CaSiO(HO)3+
0, 0.0 , 0.0 // 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;
}
void solve()
{
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 = 0.2;
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().maxstep = 200.0;
solver.get_options().maxiter_maxstep = 100;
solver.get_options().factor_descent_condition = 1e-4;
Eigen::VectorXd x(prog->total_variables());
//x.setConstant(-1);
x(0) = 1.0;
x(1) = 0.0;
x(2) = 0.0;
x(3) = 0.0;
x(4) = 0.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()
{
solve();
//solve_aluminium();
}
c4science · Help