Page MenuHomec4science
Files F73360500

equilibrium_curve.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Jul 20, 20:45

equilibrium_curve.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"
class LeachingEquilibriumCurve: public specmicp::EquilibriumCurve
{
public:
specmicp::scalar_t step = 25;
LeachingEquilibriumCurve()
{
specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
});
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 = 7e3;
specmicp::scalar_t m_c3s = mult*0.7;
specmicp::scalar_t m_c2s = mult*0.3;
specmicp::scalar_t wc = 0.5;
specmicp::scalar_t m_water = wc*1e-3*(
m_c3s*(3*56.08+60.08)
+ m_c2s*(2*56.06+60.08)
);
formulation.mass_solution = m_water;
formulation.amount_minerals = {
{"C3S", m_c3s},
{"C2S", m_c2s},
};
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_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 = 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;
solve_first_problem();
std::cout << "total concentration \t TOTaq \t TOTs \t S^t_m \t pH";
for (specmicp::index_t mineral: raw_data->range_mineral())
{
std::cout << "\t" << raw_data->get_label_mineral(mineral);
}
std::cout << std::endl;
specmicp::AdimensionalSystemSolutionExtractor sol(current_solution(), database(), solver_options().units_set);
specmicp::scalar_t tot_conc_ca = constraints().total_concentrations(id_ca);
total_steps = std::ceil(tot_conc_ca / step) + 1;
cnt = 0;
solution = specmicp::Matrix::Zero(total_steps, 2);
output();
}
void output()
{
specmicp::AdimensionalSystemSolutionExtractor sol(current_solution(), database(), solver_options().units_set);
std::cout << constraints().total_concentrations(id_ca) << "\t"
<< sol.total_aqueous_concentration(id_ca) << "\t"
<< sol.total_solid_concentration(id_ca) << "\t"
<< sol.total_saturation_minerals() << "\t"
<< sol.pH();
for (specmicp::index_t mineral: database()->range_mineral())
{
std::cout << "\t" << sol.mole_concentration_mineral(mineral);
}
std::cout << std::endl;
solution(cnt, 0) = sol.total_solid_concentration(id_ca);
solution(cnt, 1) = sol.total_aqueous_concentration(id_ca);
++cnt;
}
specmicp::index_t nb_steps() const
{
return total_steps;
}
void update_problem()
{
constraints().total_concentrations(id_ca) -= step;
constraints().total_concentrations(id_ho) += 2*step;
}
specmicp::Matrix& equilibrium_curve() {return solution;}
private:
mutable specmicp::index_t cnt;
specmicp::index_t total_steps;
mutable specmicp::Matrix solution;
specmicp::index_t id_h2o;
specmicp::index_t id_ho;
specmicp::index_t id_ca;
};
int main()
{
LeachingEquilibriumCurve solver;
for (int i=1; i<solver.nb_steps(); ++i)
{
solver.run_step();
}
std::cout << solver.equilibrium_curve() << std::endl;
}

Event Timeline

c4science · Help