Homec4science
Diffusion SpecMiCP / ReactMiCP (gasphase) Browse

specmicp/tests/specmicp/rate_nicoleau.cpp
e4061f7461edgasphase

Last Change
Raw File

rate_nicoleau.cpp
View Options

#include <iostream>
#include <fstream>
#include <vector>
#include <array>
#include <stdexcept>
#include <specmicp/reaction_path.hpp>
// Compute the supersaturation of C3S to obtain the kinetic rate curve r=r(SI)
// Data from Nicoleau et al. (2013)
//
// Nicoleau, L., Nonat, A., Perrey, D.:
// The di- and tricalcium silicate dissolutions
// Cement and Concrete Research 47(0), 14–30, 2013
const int nb_cells = 9;
const double Mc3s = 228.3;
using input_data_t = std::array<double, nb_cells>;
void read_csv_file(const std::string& filepath, std::vector<input_data_t>& vector_data)
{
std::ifstream datafile(filepath);
if (not datafile.is_open())
{
throw std::invalid_argument("Input file cannot be opened : " + filepath);
}
std::string buffer;
getline(datafile, buffer); // first line is headers
while (datafile.good())
{
input_data_t data;
for (int i= 0; i<nb_cells; ++i)
{
getline(datafile, buffer, ',');
//std::cout << buffer << std::endl;
try {
data[i] = std::stod(buffer);
}
catch (const std::invalid_argument& e)
{
throw std::invalid_argument("std : "+buffer);
}
}
vector_data.push_back(data);
}
}
double get_sursaturation(const input_data_t& input)
{
double m_nacl = 1e-3*input[2];
double m_caoh2 = 1e-3*input[3];
double m_naoh = 1e-3*input[4];
double m_cacl2 = 1e-3*input[5];
double m_c3s = input[1]/100/input[0]/Mc3s + 1e-3*input[6];
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_water = 1.0;
model->amount_aqueous = {
{"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water, 0)},
{"Na[+]", specmicp::reaction_amount_t(m_naoh+m_nacl, 0)},
{"Cl[-]", specmicp::reaction_amount_t(2*m_cacl2, 0)},
{"Ca[2+]", specmicp::reaction_amount_t(m_cacl2, 0)},
{"HO[-]", specmicp::reaction_amount_t(m_naoh, 0)}
};
model->amount_minerals = {
{"C3S", specmicp::reaction_amount_t(m_c3s, 0)},
{"Portlandite", specmicp::reaction_amount_t(m_caoh2, 0)},
};
specmicp::ReactionPathDriver driver(model, data);
driver.dissolve_to_components();
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.);
specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x);
if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized)
{
std::cout << "Error : problem not solved : return code " << (int) perf.return_code << std::endl;
}
return driver.get_current_solution().logIAP_kinetic("C3S");
}
int main()
{
std::vector<input_data_t> input;
//read_csv_file("data_test/data_nicoleau_c3sm.csv", input);
//read_csv_file("data_test/data_nicoleau_c3st1.csv", input);
read_csv_file("data_test/data_nicoleau_c3st2.csv", input);
// for (auto it=input.begin(); it!=input.end(); ++it)
// {
// for (auto its=it->begin(); its!=it->end(); ++its)
// {
// std::cout << *its << "\t";
// }
// std::cout << std::endl;
// }
for (auto it=input.begin(); it!=input.end(); ++it)
{
std::cout << get_sursaturation(*it) << "\t" << (*it)[7] << "\t" << (*it)[8] << std::endl;
}
return EXIT_SUCCESS;
}
c4science · Help