Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/tests/reactmicp/systems/saturated_react/speciation_system.cpp
a4ba90c69915reactmicp/unsaturated

Last Change
Raw File

speciation_system.cpp
View Options

#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include <iostream>
specmicp::RawDatabasePtr leaching_database()
{
specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
});
thedatabase.keep_only_components({"H2O", "H[+]", "Ca[2+]", "Si(OH)4"});
thedatabase.remove_gas_phases();
thedatabase.swap_components(swapping);
return thedatabase.get_database();
}
specmicp::AdimensionalSystemSolution initial_leaching_pb(
specmicp::RawDatabasePtr raw_data,
specmicp::scalar_t mult,
specmicp::units::UnitsSet the_units)
{
specmicp::Formulation formulation;
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*1.0e-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);
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
constraints.charge_keeper = 1;
constraints.set_saturated_system();
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 100.0;
options.solver_options.max_iter = 100;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = true;
options.solver_options.factor_descent_condition = -1.0;
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.solver_options.non_monotone_linesearch = false;
options.system_options.non_ideality_tolerance = 1e-10;
options.units_set = the_units;
specmicp::Vector x(raw_data->nb_component+raw_data->nb_mineral);
x.setZero();
x(0) = 0.5;
x.segment(1, raw_data->nb_component-1).setConstant(-3.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);
return solver.get_raw_solution(x);
}
specmicp::AdimensionalSystemSolution initial_blank_leaching_pb(
specmicp::RawDatabasePtr raw_data,
specmicp::scalar_t mult,
specmicp::units::UnitsSet the_units
)
{
specmicp::Formulation formulation;
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.concentration_aqueous = {
{"HO[-]", 1e-9},
{"Ca[2+]", 1e-9},
{"SiO(OH)3[-]", 1e-9}
};
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
constraints.charge_keeper = 1;
constraints.set_saturated_system();
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;
options.units_set = the_units;
specmicp::Vector x(raw_data->nb_component+raw_data->nb_mineral);
x.setZero();
x(0) = 0.8;
x.segment(1, raw_data->nb_component-1).setConstant(-3.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);
return solver.get_raw_solution(x);
}
c4science · Help