Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/examples/specmicp/adimensional/thermocarbo.cpp
498ace618930reactmicp/unsaturated

Last Change
Raw File

thermocarbo.cpp
View Options

#include <iostream>
#include "utils/log.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "specmicp/adimensional/equilibrium_curve.hpp"
#include "database/database.hpp"
void fixed_fugacity_thermocarbo()
{
specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"Al[3+]","Al(OH)4[-]"},
});
thedatabase.swap_components(swapping);
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::Formulation formulation;
specmicp::scalar_t mult = 5e3;
specmicp::scalar_t m_c3s = mult*0.6;
specmicp::scalar_t m_c2s = mult*0.2;
specmicp::scalar_t m_c3a = mult*0.10;
specmicp::scalar_t m_gypsum = mult*0.10;
specmicp::scalar_t wc = 0.8;
specmicp::scalar_t m_water = wc*1e-3*(
m_c3s*(3*56.08+60.08)
+ m_c2s*(2*56.06+60.08)
+ m_c3a*(3*56.08+101.96)
+ m_gypsum*(56.08+80.06+2*18.02)
);
formulation.mass_solution = m_water;
formulation.amount_minerals = {
{"C3S", m_c3s},
{"C2S", m_c2s},
{"C3A", m_c3a},
{"Gypsum", m_gypsum}
};
formulation.extra_components_to_keep = {"HCO3[-]", };
formulation.minerals_to_keep = {
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2,am",
"Calcite",
"Al(OH)3,am",
"Monosulfoaluminate",
"Tricarboaluminate",
"Monocarboaluminate",
"Hemicarboaluminate",
"Straetlingite",
"Gypsum",
"Ettringite",
"Thaumasite"
};
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
//specmicp::index_t id_h2o = thedatabase.component_label_to_id("H2O");
specmicp::index_t id_ho = thedatabase.component_label_to_id("HO[-]");
specmicp::index_t id_hco3 = thedatabase.component_label_to_id("HCO3[-]");
specmicp::index_t id_co2g = thedatabase.gas_label_to_id("CO2(g)");
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
constraints.charge_keeper = id_ho;
constraints.add_fixed_fugacity_gas(id_co2g, id_hco3, -14);
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 50.0;
options.solver_options.max_iter = 50;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = false;
options.solver_options.factor_descent_condition = -1;
options.solver_options.factor_gradient_search_direction = 100;
options.solver_options.projection_min_variable = 1e-6;
options.solver_options.fvectol = 1e-8;
options.solver_options.steptol = 1e-10;
options.solver_options.non_monotone_linesearch = false;
options.system_options.non_ideality_tolerance = 1e-14;
specmicp::Vector x(raw_data->nb_component()+raw_data->nb_mineral());
x.setZero();
x(0) = 0.5;
x.segment(2, raw_data->nb_component()).setConstant(-6.0);
for (specmicp::index_t i: raw_data->range_component())
{
std::cout << raw_data->components.get_label(i) << std::endl;
}
specmicp::uindex_t tot_iter = 0;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
specmicp_assert((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
//std::cout << perf.nb_iterations << std::endl;
tot_iter += perf.nb_iterations;
std::cout << "log fugacity" << "\t" << "pH";
for (specmicp::index_t mineral: raw_data->range_mineral())
{
std::cout << "\t" << raw_data->minerals.get_label(mineral);
}
std::cout << std::endl;
specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor sol(solution, raw_data, options.units_set);
std::cout << -14 << "\t" << sol.pH();
for (specmicp::index_t mineral: raw_data->range_mineral())
{
std::cout << "\t" << sol.mole_concentration_mineral(mineral);
}
std::cout << std::endl;
for (double logf=-13.9; logf<-4.1; logf+=0.05)
{
constraints.fixed_fugacity_cs[0].log_value = logf;
solver = specmicp::AdimensionalSystemSolver (raw_data, constraints, solution, options);
//std::cout << solver.get_options().solver_options.factor_descent_condition << std::endl;
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
specmicp_assert((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
std::cout << perf.nb_iterations << std::endl;
tot_iter += perf.nb_iterations;
solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor sol(solution, raw_data, options.units_set);
std::cout << logf << "\t" << sol.pH();
for (specmicp::index_t mineral: raw_data->range_mineral())
{
std::cout << "\t" << sol.mole_concentration_mineral(mineral);
}
std::cout << std::endl;
//specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
//std::cout << solution.main_variables(0) << std::endl;
//std::cout << 14+solution.main_variables(1) << std::endl;
}
std::cout << tot_iter << std::endl;
}
void reaction_path_thermocarbo()
{
specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"Al[3+]","Al(OH)4[-]"},
});
thedatabase.swap_components(swapping);
thedatabase.remove_gas_phases();
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::Formulation formulation;
specmicp::scalar_t mult = 5e3;
specmicp::scalar_t m_c3s = mult*0.6;
specmicp::scalar_t m_c2s = mult*0.2;
specmicp::scalar_t m_c3a = mult*0.10;
specmicp::scalar_t m_gypsum = mult*0.10;
specmicp::scalar_t wc = 0.8;
specmicp::scalar_t m_water = wc*1e-3*(
m_c3s*(3*56.08+60.08)
+ m_c2s*(2*56.06+60.08)
+ m_c3a*(3*56.08+101.96)
+ m_gypsum*(56.08+80.06+2*18.02)
);
formulation.mass_solution = m_water;
formulation.amount_minerals = {
{"C3S", m_c3s},
{"C2S", m_c2s},
{"C3A", m_c3a},
{"Gypsum", m_gypsum}
};
formulation.extra_components_to_keep = {"HCO3[-]", };
formulation.minerals_to_keep = {
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2,am",
"Calcite",
"Al(OH)3,am",
"Monosulfoaluminate",
"Tricarboaluminate",
"Monocarboaluminate",
"Hemicarboaluminate",
"Straetlingite",
"Gypsum",
"Ettringite",
"Thaumasite"
};
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
specmicp::index_t id_h2o = thedatabase.component_label_to_id("H2O");
specmicp::index_t id_ho = thedatabase.component_label_to_id("HO[-]");
specmicp::index_t id_hco3 = thedatabase.component_label_to_id("HCO3[-]");
//specmicp::index_t id_co2g = thedatabase.gas_label_to_id("CO2(g)");
specmicp::index_t id_ca = thedatabase.component_label_to_id("Ca[2+]");
std::cout << total_concentrations << std::endl;
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
//constraints.charge_keeper = id_ho;
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 10.0;
options.solver_options.max_iter = 100;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = false;
options.solver_options.factor_descent_condition = -1;
options.solver_options.factor_gradient_search_direction = 100;
options.solver_options.projection_min_variable = 1e-9;
options.solver_options.fvectol = 1e-6;
options.solver_options.steptol = 1e-14;
options.system_options.non_ideality_tolerance = 1e-10;
specmicp::Vector x(raw_data->nb_component()+raw_data->nb_mineral());
x.setZero();
x(0) = 0.8;
x.segment(2, raw_data->nb_component()).setConstant(-6.0);
for (specmicp::index_t i: raw_data->range_component())
{
std::cout << raw_data->components.get_label(i) << std::endl;
}
specmicp::uindex_t tot_iter = 0;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
specmicp_assert((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
//std::cout << perf.nb_iterations << std::endl;
tot_iter += perf.nb_iterations;
std::cout << "Buffering \t" << "pH";
for (specmicp::index_t mineral: raw_data->range_mineral())
{
std::cout << "\t" << raw_data->minerals.get_label(mineral);
}
std::cout << std::endl;
specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor sol(solution, raw_data, options.units_set);
std::cout << 0 << "\t" << sol.pH();
for (specmicp::index_t mineral: raw_data->range_mineral())
{
std::cout << "\t" << sol.mole_concentration_mineral(mineral);
}
std::cout << std::endl;
for (double c_hco3=100; c_hco3<=constraints.total_concentrations(id_ca); c_hco3+=100)
{
constraints.total_concentrations(id_hco3) = c_hco3;
constraints.total_concentrations(id_ho) -= 100;
constraints.total_concentrations(id_h2o) += 100;
solver = specmicp::AdimensionalSystemSolver (raw_data, constraints, solution, options);
//std::cout << solver.get_options().solver_options.factor_descent_condition << std::endl;
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
specmicp_assert((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
//std::cout << perf.nb_iterations << std::endl;
tot_iter += perf.nb_iterations;
solution = solver.get_raw_solution(x);
specmicp::AdimensionalSystemSolutionExtractor sol(solution, raw_data, options.units_set);
std::cout << c_hco3/constraints.total_concentrations(id_ca) << "\t" << sol.pH();
for (specmicp::index_t mineral: raw_data->range_mineral())
{
std::cout << "\t" << sol.mole_concentration_mineral(mineral);
}
std::cout << std::endl;
//specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
//std::cout << solution.main_variables(0) << std::endl;
//std::cout << 14+solution.main_variables(1) << std::endl;
}
std::cout << tot_iter << std::endl;
}
class AutomaticReactionPath: public specmicp::EquilibriumCurve
{
public:
AutomaticReactionPath()
{
specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"Al[3+]","Al(OH)4[-]"},
});
thedatabase.swap_components(swapping);
thedatabase.remove_gas_phases();
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
set_database(raw_data);
specmicp::Formulation formulation;
specmicp::scalar_t mult = 5e3;
specmicp::scalar_t m_c3s = mult*0.6;
specmicp::scalar_t m_c2s = mult*0.2;
specmicp::scalar_t m_c3a = mult*0.10;
specmicp::scalar_t m_gypsum = mult*0.10;
specmicp::scalar_t wc = 0.8;
specmicp::scalar_t m_water = wc*1e-3*(
m_c3s*(3*56.08+60.08)
+ m_c2s*(2*56.06+60.08)
+ m_c3a*(3*56.08+101.96)
+ m_gypsum*(56.08+80.06+2*18.02)
);
formulation.mass_solution = m_water;
formulation.amount_minerals = {
{"C3S", m_c3s},
{"C2S", m_c2s},
{"C3A", m_c3a},
{"Gypsum", m_gypsum}
};
formulation.extra_components_to_keep = {"HCO3[-]", };
formulation.minerals_to_keep = {
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2,am",
"Calcite",
"Al(OH)3,am",
"Monosulfoaluminate",
"Tricarboaluminate",
"Monocarboaluminate",
"Hemicarboaluminate",
"Straetlingite",
"Gypsum",
"Ettringite",
"Thaumasite"
};
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
id_h2o = thedatabase.component_label_to_id("H2O");
id_ho = thedatabase.component_label_to_id("HO[-]");
id_hco3 = thedatabase.component_label_to_id("HCO3[-]");
id_co2g = thedatabase.gas_label_to_id("CO2(g)");
id_ca = thedatabase.component_label_to_id("Ca[2+]");
constraints() = specmicp::AdimensionalSystemConstraints(total_concentrations);
//constraints().charge_keeper = id_ho;
specmicp::AdimensionalSystemSolverOptions& options = solver_options();
options.solver_options.maxstep = 10.0;
options.solver_options.max_iter = 200;
options.solver_options.maxiter_maxstep = 200;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = false;
options.solver_options.non_monotone_linesearch = false;
options.solver_options.factor_descent_condition = -1;
options.solver_options.factor_gradient_search_direction = 100;
options.solver_options.projection_min_variable = 1e-9;
options.solver_options.fvectol = 1e-8;
options.solver_options.steptol = 1e-10;
options.system_options.non_ideality_tolerance = 1e-10;
solve_first_problem();
std::cout << "Buffering \t" << "pH";
for (specmicp::index_t mineral: raw_data->range_mineral())
{
std::cout << "\t" << raw_data->minerals.get_label(mineral);
}
std::cout << std::endl;
output();
}
void output()
{
specmicp::AdimensionalSystemSolutionExtractor sol(current_solution(), database(), solver_options().units_set);
std::cout << constraints().total_concentrations(id_hco3)/constraints().total_concentrations(id_ca) << "\t" << sol.pH();
for (specmicp::index_t mineral: database()->range_mineral())
{
std::cout << "\t" << sol.mole_concentration_mineral(mineral);
}
std::cout << std::endl;
}
void update_problem()
{
constraints().total_concentrations(id_hco3) += 100;
constraints().total_concentrations(id_ho) -= 100;
constraints().total_concentrations(id_h2o) += 100;
}
private:
specmicp::index_t id_h2o;
specmicp::index_t id_ho;
specmicp::index_t id_hco3;
specmicp::index_t id_co2g;
specmicp::index_t id_ca;
};
int main()
{
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
//fixed_fugacity_thermocarbo();
//reaction_path_thermocarbo();
AutomaticReactionPath test_automatic;
for (int i=0; i<130; ++i)
{
test_automatic.run_step();
}
return EXIT_SUCCESS;
}
c4science · Help