Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/examples/specmicp/adimensional/momas_thermo.cpp
8c94dee2ed89reactivetransport

Last Change
Raw File

momas_thermo.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 "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/adimensional_kinetics/kinetic_variables.hpp"
#include "specmicp/adimensional_kinetics/kinetic_model.hpp"
#include "specmicp/adimensional_kinetics/kinetic_system_euler_solver.hpp"
#include "specmicp/adimensional_kinetics/kinetic_system_solver.hpp"
specmicp::RawDatabasePtr get_momas_db()
{
specmicp::database::Database thedatabase("../data/momas_benchmark.js");
thedatabase.remove_components({thedatabase.component_label_to_id("X5")});
return thedatabase.get_database();
}
class AutomaticReactionPath: public specmicp::EquilibriumCurve
{
public:
AutomaticReactionPath()
{
specmicp::RawDatabasePtr raw_data = get_momas_db();
set_database(raw_data);
std::cout << raw_data->nu_aqueous << std::endl;
std::cout << raw_data->logk_aqueous << std::endl;
std::cout << raw_data->nu_mineral << std::endl;
std::cout << raw_data->logk_mineral << std::endl;
specmicp::Vector total_concentrations(6);
total_concentrations << 1, 0.0, -3.0, 0.0, 1.0, 1.0;
constraints() = specmicp::AdimensionalSystemConstraints(total_concentrations);
constraints().disable_conservation_water();
specmicp::AdimensionalSystemSolverOptions& options = solver_options();
options.solver_options.maxstep = 1.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;
options.system_options.non_ideality = false;
//solve_first_problem();
specmicp::Vector variables(raw_data->nb_component+raw_data->nb_mineral);
variables << 0.8, -5, -1, -5, -1, -3, 0.0;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints(), solver_options());
specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, false);
std::cout << (int) current_perf.return_code << std::endl;
solution_vector() = variables;
initialize_solution(solver.get_raw_solution(variables));
std::cout << "Buffering \t";
for (specmicp::index_t mineral: raw_data->range_mineral())
{
std::cout << "\t" << raw_data->labels_minerals[mineral];
}
for (specmicp::index_t sorbed: raw_data->range_sorbed())
{
std::cout << "\t" << raw_data->labels_sorbed[sorbed];
}
std::cout << std::endl;
output();
}
void output()
{
specmicp::AdimensionalSystemSolutionExtractor sol(current_solution(), database(), solver_options().units_set);
std::cout << constraints().total_concentrations(1);
for (specmicp::index_t mineral: database()->range_mineral())
{
std::cout << "\t" << sol.mole_concentration_mineral(mineral);
}
for (specmicp::index_t sorbed: database()->range_sorbed())
{
std::cout << "\t" << sol.sorbed_species_molalities(sorbed);
}
std::cout << std::endl;
}
void update_problem()
{
constraints().total_concentrations(1) += 0.1;
constraints().total_concentrations(2) += 0.1;
constraints().total_concentrations(3) += 0.1;
}
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;
};
class MomasKinetics: public specmicp::kinetics::AdimKineticModel
{
public:
MomasKinetics():
specmicp::kinetics::AdimKineticModel({0}),
m_id_cc(0),
m_id_c3(2),
m_id_x4(4)
{}
//! \brief Compute the kinetic rates and store them in dydt
void compute_rate(
specmicp::scalar_t t,
const specmicp::Vector& y,
specmicp::kinetics::AdimKineticVariables& variables,
specmicp::Vector& dydt
)
{
dydt.resizeLike(y);
specmicp::AdimensionalSystemSolution& solution = variables.equilibrium_solution();
specmicp::scalar_t factor = std::pow(1e3*0.5*solution.secondary_molalities(m_id_c3),3)/
std::pow(1e3*0.5*pow10(solution.main_variables(m_id_x4)), 2);
factor = 0.2*factor-1.0;
specmicp::scalar_t k_constant = 10.0;
if (factor >= 0)
k_constant = 1e-2;
dydt(0) = k_constant*factor;
}
private:
specmicp::index_t m_id_cc;
specmicp::index_t m_id_c3;
specmicp::index_t m_id_x4;
};
void solve_kinetics_problem()
{
specmicp::RawDatabasePtr raw_data = get_momas_db();
specmicp::Vector total_concentrations(6);
total_concentrations << 1, 0.6, -2.4, 0.6, 1.0, 1.0;
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
constraints.disable_conservation_water();
specmicp::AdimensionalSystemSolverOptions 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 = true;
options.solver_options.non_monotone_linesearch = false;
options.solver_options.factor_descent_condition = -1;
options.solver_options.factor_gradient_search_direction = 50;
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;
options.system_options.non_ideality = false;
specmicp::Vector equil_variables(raw_data->nb_component+raw_data->nb_mineral);
equil_variables << 0.8, -5, -1, -5, -1, -3, 0.0;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(equil_variables, false);
std::cout << (int) current_perf.return_code << std::endl;
specmicp::AdimensionalSystemSolution equil_solution = solver.get_raw_solution(equil_variables);
std::shared_ptr<specmicp::kinetics::AdimKineticModel> model = std::make_shared<MomasKinetics>();
specmicp::Vector mineral_concentrations(1);
mineral_concentrations << 5.0;
specmicp::kinetics::AdimKineticSystemEulerSolver kinetic_solver(
//specmicp::kinetics::AdimKineticSystemSolver kinetic_solver(
model,
total_concentrations,
mineral_concentrations,
constraints,
equil_solution,
raw_data);
options.solver_options.use_scaling = false;
kinetic_solver.get_options().speciation_options = options;
kinetic_solver.solve(0.005, 1.0);
std::cout << "Final amount : " << kinetic_solver.variables().concentration_mineral(0) << std::endl;
}
int main()
{
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
solve_kinetics_problem();
//fixed_fugacity_thermocarbo();
//reaction_path_thermocarbo();
// AutomaticReactionPath test_automatic;
// for (int i=0; i<20; ++i)
// {
// test_automatic.run_step();
// }
return EXIT_SUCCESS;
}
c4science · Help