Homec4science
Diffusion SpecMiCP / ReactMiCP (dev) Browse

specmicp/tests/reactmicp/systems/diffusion/utils.cpp
afebadac83dedev

Last Change
Raw File

utils.cpp
View Options

#include "utils.hpp"
EquilibriumState blank_composition(std::shared_ptr<specmicp::database::DataContainer> database)
{
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
double m_c3s = 1e-7;
double m_c2s = 1e-7;
double m_water = 1;
double c_hco3 = 1e-4;
model->amount_aqueous = {
{"H2O", specmicp::reaction_amount_t(m_water/database->molar_mass_basis_si(0)+c_hco3, 0)},
{"HCO3[-]", specmicp::reaction_amount_t(c_hco3, 0)},
{"HO[-]", specmicp::reaction_amount_t(-c_hco3, 0)},
};
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 x;
specmicp::ReactionPathDriver driver(model, database, x);
//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;
specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x);
if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized)
throw std::runtime_error("Error - system cannot be solved");
return driver.get_current_solution();
}
EquilibriumState sample_carbo_composition(std::shared_ptr<specmicp::database::DataContainer> database)
{
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
const double mult = 0.005;
const double m_c3s = mult*0.7;
const double m_c2s = mult*0.3;
const double wc = 1.0;
const double m_water = wc*((3*56.08+60.08)*m_c3s + (2*56.08+60.08)*m_c2s)*1e-3;
const double delta_h2co3 = 0.005;
model->amount_aqueous = {
{"H2O", specmicp::reaction_amount_t(m_water/database->molar_mass_basis_si(0), delta_h2co3)},
{"HCO3[-]", specmicp::reaction_amount_t(delta_h2co3, 0.0)},
{"HO[-]", specmicp::reaction_amount_t(-delta_h2co3, 0.0)},
};
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 x;
specmicp::ReactionPathDriver driver(model, database, x);
//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;
specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x);
if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized)
throw std::runtime_error("Error - system cannot be solved");
return driver.get_current_solution();
}
//! \brief Return a simple fully processed database for testing purposes
//!
//! Correspond to the carbonation problem
std::shared_ptr<specmicp::database::DataContainer> get_test_carbo_database()
{
specmicp::database::Database database("data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"}
});
database.swap_components(swapping);
std::vector<std::string> to_keep = {"H2O", "HO[-]", "Ca[2+]", "SiO(OH)3[-]", "HCO3[-]"};
database.keep_only_components(to_keep);
return database.get_database();
}
c4science · Help