Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F68274663
momas_thermo.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Jun 26, 16:35
Size
7 KB
Mime Type
text/x-c++
Expires
Fri, Jun 28, 16:35 (2 d)
Engine
blob
Format
Raw Data
Handle
18554173
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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;
}
Event Timeline
Log In to Comment