Homec4science
Diffusion SpecMiCP / ReactMiCP (bin/unsaturated) Browse

specmicp/examples/reactmicp/saturated_react/react_leaching.cpp
2553b8ca7134bin/unsaturated

Last Change
Raw File

react_leaching.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 "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp"
#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
#include "reactmicp/solver/reactive_transport_solver.hpp"
#include "reactmicp/systems/saturated_react/variables.hpp"
#include "reactmicp/systems/saturated_react/transport_stagger.hpp"
#include "reactmicp/systems/saturated_react/equilibrium_stagger.hpp"
#include "reactmicp/systems/saturated_react/init_variables.hpp"
#include "dfpm/meshes/axisymmetric_uniform_mesh1d.hpp"
#include "dfpm/meshes/axisymmetric_mesh1d.hpp"
#include "database/database.hpp"
#include "reactmicp/solver/timestepper.hpp"
#include "utils/timer.hpp"
#include "utils/log.hpp"
#include <iostream>
#include <fstream>
#include <Eigen/LU>
#include "utils/io/reactive_transport.hpp"
using namespace specmicp;
// Initial conditions
// ==================
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[-]"},
{"Al[3+]","Al(OH)4[-]"}
});
thedatabase.remove_gas_phases();
thedatabase.swap_components(swapping);
return thedatabase.get_database();
}
const double M_CaO = 56.08;
const double M_SiO2 = 60.09;
const double M_Al2O3 = 101.96;
const double M_SO3 = 80.06;
void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species)
{
Eigen::MatrixXd compo_in_oxyde(4, 4);
Vector molar_mass(4);
molar_mass << 1.0/M_CaO, 1.0/M_SiO2, 1.0/M_Al2O3, 1.0/M_SO3;
compo_oxyde = compo_oxyde.cwiseProduct(molar_mass);
//C3S C2S C3A Gypsum
compo_in_oxyde << 3, 2, 3, 1, // CaO
1, 1, 0, 0, // SiO2
0, 0, 1, 0, // Al2O3
0, 0, 0, 1; // SO3
Eigen::FullPivLU<Eigen::MatrixXd> solver(compo_in_oxyde);
compo_species = solver.solve(compo_oxyde);
}
AdimensionalSystemSolution initial_leaching_pb(
RawDatabasePtr raw_data,
scalar_t mult,
units::UnitsSet the_units)
{
Formulation formulation;
Vector oxyde_compo(4);
oxyde_compo << 62.9, 20.6, 5.8, 3.1;
Vector species_compo;
compo_from_oxyde(oxyde_compo, species_compo);
std::cout << species_compo << std::endl;
scalar_t m_c3s = mult*species_compo(0); //0.6;
scalar_t m_c2s = mult*species_compo(1); //0.2;
scalar_t m_c3a = mult*species_compo(2); //0.1;
scalar_t m_gypsum = mult*species_compo(3); //0.1;
scalar_t wc = 0.5;
scalar_t m_water = wc*1.0e-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}
};
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 = 200;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.projection_min_variable = 1e-9;
options.solver_options.set_tolerance(1e-10, 1e-14);
options.solver_options.disable_crashing();
options.solver_options.enable_scaling();
options.solver_options.disable_descent_direction();
options.solver_options.enable_non_monotone_linesearch();
options.system_options.non_ideality_tolerance = 1e-10;
options.units_set = the_units;
specmicp::Vector x;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
solver.initialise_variables(x, 0.8, -3.0);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
std::cout << "Initial return code : " << (int) perf.return_code << std::endl;
if ((int) perf.return_code <= (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet)
{
throw std::runtime_error("Initial solution has not converged");
}
return solver.get_raw_solution(x);
}
AdimensionalSystemSolution initial_blank_leaching_pb(
RawDatabasePtr raw_data,
scalar_t mult,
units::UnitsSet the_units)
{
Formulation formulation;
Vector oxyde_compo(4);
oxyde_compo << 62.9, 20.6, 5.8, 3.1;
Vector species_compo;
compo_from_oxyde(oxyde_compo, species_compo);
std::cout << species_compo << std::endl;
scalar_t m_c3s = 0.6;
scalar_t m_c2s = 0.2;
scalar_t m_c3a = 0.1;
scalar_t m_gypsum = 0.1;
scalar_t wc = 0.8;
scalar_t m_water = wc*1.0e-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 = {
{"SiO2,am", mult},
};
// formulation.concentration_aqueous = {
// {"Ca[2+]", 1e-7},
// {"Al(OH)4[-]", 1e-7},
// {"SO4[2-]", 1e-7}
// };
formulation.extra_components_to_keep = {"Ca[2+]", "SO4[2-]", "Al(OH)4[-]"};
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 = 200;
options.solver_options.maxiter_maxstep = 200;
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-4;
options.solver_options.steptol = 1e-14;
options.solver_options.non_monotone_linesearch = true;
options.system_options.non_ideality_tolerance = 1e-10;
options.units_set = the_units;
Vector x;
AdimensionalSystemSolver solver(raw_data, constraints, options);
solver.initialise_variables(x, 0.75, -4.0);
micpsolver::MiCPPerformance perf = solver.solve(x, false);
std::cout << "Initial return code : " << (int) perf.return_code << std::endl;
if ((int) perf.return_code <= (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet)
{
throw std::runtime_error("Initial solution has not converged");
}
AdimensionalSystemSolution solution = solver.get_raw_solution(x);
AdimensionalSystemSolutionModificator extractor(solution, raw_data, the_units);
extractor.remove_solids();
Vector totconc(raw_data->nb_component());
totconc(0) = extractor.density_water()*extractor.total_aqueous_concentration(0);
for (index_t component: raw_data->range_aqueous_component())
{
totconc(component) = 0.01*extractor.density_water()*extractor.total_aqueous_concentration(component);
}
constraints.total_concentrations = totconc;
solver = AdimensionalSystemSolver(raw_data, constraints, options);
perf = solver.solve(x, false);
specmicp_assert((int) perf.return_code > (int) micpsolver::MiCPSolverReturnCode::NotConvergedYet);
return solver.get_raw_solution(x);
}
// Upscaling
// =========
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::solver;
using namespace specmicp::reactmicp::systems::satdiff;
class LeachingUspcaling: public solver::UpscalingStaggerBase
{
public:
LeachingUspcaling(RawDatabasePtr the_database,
units::UnitsSet the_units):
m_data(the_database),
m_units_set(the_units)
{
}
//! \brief Initialize the stagger at the beginning of the computation
void initialize(VariablesBasePtr var)
{
SaturatedVariablesPtr true_var = cast_var_from_base(var);
database::Database db_handler(m_data);
m_dt = 1;
m_id_cshj = db_handler.mineral_label_to_id("CSH,j");
m_initial_sat_cshj = true_var->equilibrium_solution(1).main_variables(m_data->nb_component()+m_id_cshj);
//m_initial_sat_cshj = 0.36;
std::cout << m_initial_sat_cshj << std::endl;
true_var->upscaling_variables().setZero();
for (index_t node=0; node<true_var->nb_nodes(); ++node)
{
upscaling_one_node(node, true_var);
}
}
//! \brief Initialize the stagger at the beginning of an iteration
void initialize_timestep(scalar_t dt, VariablesBasePtr var)
{
SaturatedVariablesPtr true_var = cast_var_from_base(var);
m_dt = dt;
for (index_t node=1; node<true_var->nb_nodes(); ++node)
{
upscaling_one_node(node, true_var);
true_var->vel_porosity(node) = 0.0;
}
}
//! \brief Solve the equation for the timestep
StaggerReturnCode restart_timestep(VariablesBasePtr var)
{
SaturatedVariablesPtr true_var = cast_var_from_base(var);
for (index_t node=1; node<true_var->nb_nodes(); ++node)
{
upscaling_one_node(node, true_var);
}
return StaggerReturnCode::ResidualMinimized;
}
scalar_t diffusion_coefficient(scalar_t porosity)
{
const scalar_t factor = 1e4;
scalar_t poro = factor*std::exp(9.85*porosity-29.08);
// const scalar_t de_0 = 2.76e-12;
// const scalar_t porosity_r = 0.02;
// const scalar_t exponent = 3.32;
//if (porosity > 0.7)
// return factor*4.0*1e-10;
// else
// {
// const scalar_t ratio = ((porosity - porosity_r)/(0.25 - porosity_r));
// poro = factor*de_0*std::pow(ratio, exponent);
// }
return std::min(poro,factor*1e-10);
}
void upscaling_one_node(index_t node, SaturatedVariablesPtr true_var)
{
AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node),
m_data, m_units_set);
scalar_t porosity = extractor.porosity();
;//- true_var->equilibrium_solution(node).main_variables(m_data->nb_component)/m_initial_sat_cshj*0.05;
//std::cout << "porosity : " << porosity
// << " - saturation water" << true_var->equilibrium_solution(node).main_variables(0) << std::endl;
true_var->vel_porosity(node) += (porosity - true_var->porosity(node))/m_dt;
true_var->porosity(node) = porosity;
true_var->diffusion_coefficient(node) = diffusion_coefficient(porosity);
}
private:
RawDatabasePtr m_data;
units::UnitsSet m_units_set;
index_t m_id_cshj;
scalar_t m_initial_sat_cshj;
scalar_t m_dt;
};
//
void set_specmicp_options(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 = true;
options.solver_options.factor_descent_condition = -1;
options.solver_options.factor_gradient_search_direction = 100;
options.solver_options.projection_min_variable = 1e-16;
options.solver_options.fvectol = 1e-10;
options.solver_options.steptol = 1e-14;
options.solver_options.non_monotone_linesearch = true;
options.system_options.non_ideality_tolerance = 1e-12;
options.system_options.non_ideality_max_iter = 100;
options.solver_options.threshold_cycling_linesearch = 1e-5;
}
void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options)
{
transport_options.maximum_iterations = 450;
transport_options.quasi_newton = 3;
transport_options.residuals_tolerance = 1e-6;
transport_options.step_tolerance = 1e-14;
transport_options.alpha = 1;
transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
transport_options.sparse_solver = SparseSolver::SparseLU;
//transport_options.sparse_solver = SparseSolver::GMRES;
transport_options.threshold_stationary_point = 1e-2;
transport_options.absolute_tolerance = 1e-16;
}
void set_reactive_solver_options(ReactiveTransportOptions& solver_options)
{
solver_options.maximum_iterations = 100;
solver_options.residuals_tolerance = 1e-2;
solver_options.good_enough_tolerance = 1.0;
solver_options.absolute_residuals_tolerance = 1e-16;
solver_options.implicit_upscaling = true;
}
mesh::Mesh1DPtr get_mesh(scalar_t dx, index_t nb_nodes, index_t additional_nodes)
{
scalar_t radius = 3.5; //cm
scalar_t height = 8.0; //cm
radius = radius + additional_nodes*dx;
return mesh::uniform_axisymmetric_mesh1d(nb_nodes+additional_nodes, radius, dx, height);
}
mesh::Mesh1DPtr get_mesh()
{
Vector radius(66);
const scalar_t height = 8.0;
radius <<
3.5005,
3.5000,
3.4995,
3.4990,
3.4980,
3.4970,
3.4960,
3.4950,
3.4925,
3.4900,
3.4875,
3.4850,
3.4825,
3.4800,
3.4775,
3.4750,
3.4725,
3.4700,
3.4675,
3.4650,
3.4625,
3.4600,
3.4575,
3.4550,
3.4525,
3.4500,
3.4475,
3.445,
3.4425,
3.440,
3.4375,
3.435,
3.4325,
3.430,
3.4275,
3.425,
3.4225,
3.420,
3.4175,
3.415,
3.4125,
3.410,
3.4075,
3.405,
3.4025,
3.400,
3.3975,
3.395,
3.3925,
3.390,
3.3875,
3.385,
3.3825,
3.38,
3.3775,
3.3750,
3.3725,
3.3700,
3.3675,
3.365,
3.3625,
3.36,
3.3575,
3.355,
3.325,
3.3
;
return mesh::axisymmetric_mesh1d(radius, height);
}
void output_mesh(mesh::Mesh1DPtr the_mesh)
{
std::ofstream output_mesh;
output_mesh.open("out_leaching_mesh.dat");
output_mesh << "Node\tPosition" << std::endl;
for (index_t node: the_mesh->range_nodes())
{
output_mesh << node << "\t" << the_mesh->get_position(node) << std::endl;
}
output_mesh.close();
}
void print_minerals_profile(
RawDatabasePtr the_database,
SaturatedVariablesPtr variables,
mesh::Mesh1DPtr the_mesh,
const std::string& filepath
)
{
std::ofstream ofile(filepath);
ofile << "Radius";
for (index_t mineral: the_database->range_mineral())
{
ofile << "\t" << the_database->get_label_mineral(mineral);
}
ofile << "\t" << "Porosity";
ofile << "\t" << "pH";
ofile << std::endl;
for (index_t node: the_mesh->range_nodes())
{
ofile << the_mesh->get_position(node);
AdimensionalSystemSolutionExtractor extractor(variables->equilibrium_solution(node), the_database, units::UnitsSet());
for (index_t mineral: the_database->range_mineral())
{
ofile << "\t" << extractor.volume_fraction_mineral(mineral);
}
ofile << "\t" << variables->porosity(node);
ofile << "\t" << extractor.pH();
ofile << std::endl;
}
ofile.close();
}
// Main
// ====
int main()
{
Eigen::initParallel();
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
RawDatabasePtr raw_data = leaching_database();
units::UnitsSet the_units;
the_units.length = units::LengthUnit::centimeter;
AdimensionalSystemSolution initial = initial_leaching_pb(raw_data, 0.015, the_units);
AdimensionalSystemSolutionModificator extractor(initial, raw_data, the_units);
index_t id_ca = raw_data->get_id_component("Ca[2+]");
index_t id_si = raw_data->get_id_component("SiO(OH)3[-]");
std::cout << extractor.total_solid_concentration(id_ca) << " - " << extractor.porosity() << std::endl;
extractor.scale_total_concentration(id_ca, 0.0145);
std::cout << extractor.total_solid_concentration(id_ca) << " - " << extractor.porosity() << std::endl;
std::cout << initial.main_variables << std::endl;
AdimensionalSystemSolution blank = initial_blank_leaching_pb(raw_data, 0.005, the_units);
std::cout << blank.main_variables << std::endl;
scalar_t dx = 0.0005;
index_t nb_nodes = 150;
index_t additional_nodes = 1;
//mesh::Mesh1DPtr the_mesh = get_mesh(dx, nb_nodes, additional_nodes);
mesh::Mesh1DPtr the_mesh = get_mesh();
nb_nodes = the_mesh->nb_nodes();
std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition = {initial, blank};
std::vector<int> index_initial_state(nb_nodes, 0);
for (int i=0; i< additional_nodes; ++i)
index_initial_state[i] = 1;
std::vector<specmicp::index_t> list_fixed_nodes = {0, };
systems::satdiff::SaturatedVariablesPtr variables =
systems::satdiff::init_variables(the_mesh, raw_data, the_units,
list_fixed_nodes,
list_initial_composition,
index_initial_state);
specmicp::AdimensionalSystemConstraints constraints;
constraints.charge_keeper = 1;
constraints.set_saturated_system();
constraints.inert_volume_fraction = 0.0;
AdimensionalSystemSolverOptions options;
set_specmicp_options(options);
options.units_set = the_units;
ChemistryStaggerPtr equilibrium_stagger =
std::make_shared<reactmicp::systems::satdiff::EquilibriumStagger>(nb_nodes, constraints, options);
std::shared_ptr<reactmicp::systems::satdiff::SaturatedTransportStagger> transport_stagger =
std::make_shared<reactmicp::systems::satdiff::SaturatedTransportStagger>(variables, list_fixed_nodes);
dfpmsolver::ParabolicDriverOptions& transport_options = transport_stagger->options_solver();
set_transport_options(transport_options);
UpscalingStaggerPtr upscaling_stagger =
std::make_shared<LeachingUspcaling>(raw_data, the_units);
ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger);
transport_stagger->initialize(variables);
equilibrium_stagger->initialize(variables);
upscaling_stagger->initialize(variables);
set_reactive_solver_options(solver.get_options());
std::ofstream output_iter("out_leaching_iter.dat");
std::ofstream output_c_ca, output_s_ca, output_s_si;
std::ofstream output_poro;
std::ofstream output_loss_ca;
output_c_ca.open("out_c_ca.dat");
output_s_ca.open("out_s_ca.dat");
output_s_si.open("out_s_si.dat");
output_poro.open("out_poro.dat");
output_loss_ca.open("out_loss_ca.dat");
scalar_t initial_ca = 0;
for (index_t node: the_mesh->range_nodes())
{
initial_ca += variables->solid_concentration(node, id_ca, variables->displacement())
* the_mesh->get_volume_cell(node);
}
scalar_t dt=2.0;
io::print_reactmicp_performance_long_header(&output_iter);
Timestepper timestepper(1.0, 500.0, 25*24*3600, 2);
int cnt =0;
while (timestepper.get_total() < timestepper.get_total_target())
{
Timer step_timer;
step_timer.start();
reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables);
step_timer.stop();
io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total() , solver.get_perfs());
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
{
dt = 1.0;
variables->reset_main_variables();
retcode = solver.solve_timestep(dt, variables);
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total(), solver.get_perfs());
}
// std::cout << "time : " << total/3600/24 << std::endl
// << "Iter : " << solver.get_perfs().nb_iterations << std::endl
// << "Residuals : " << solver.get_perfs().residuals << std::endl;
if (cnt % 100 == 0)
{
scalar_t time = timestepper.get_total()/3600.0/24.0;
scalar_t sqrt_time = std::sqrt(time);
output_c_ca << time;
output_s_ca << time;
output_s_si << time;
output_loss_ca << time << "\t" << sqrt_time;
output_poro << time;
scalar_t amount_ca = 0.0;
for (index_t node: the_mesh->range_nodes())
{
amount_ca += variables->solid_concentration(node, id_ca, variables->displacement())
* the_mesh->get_volume_cell(node);
output_c_ca << "\t" << variables->aqueous_concentration(node, id_ca, variables->displacement());
output_s_ca << "\t" << variables->solid_concentration( node, id_ca, variables->displacement());
output_s_si << "\t" << variables->solid_concentration( node, id_si, variables->displacement());
output_poro << "\t" << variables->porosity(node);
}
output_poro << std::endl;
output_loss_ca << "\t" << (initial_ca - amount_ca)/1.75929e-2 << std::endl;
output_c_ca << std::endl;
output_s_ca << std::endl;
output_s_si << std::endl;
}
++cnt;
}
print_minerals_profile(raw_data, variables, the_mesh, "out_leaching_minerals_profile_25d.dat");
timestepper.set_total_target(90*24*3600);
while (timestepper.get_total() < timestepper.get_total_target())
{
Timer step_timer;
step_timer.start();
reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables);
step_timer.stop();
io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total() , solver.get_perfs());
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
{
dt = 1.0;
variables->reset_main_variables();
retcode = solver.solve_timestep(dt, variables);
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total(), solver.get_perfs());
}
// std::cout << "time : " << total/3600/24 << std::endl
// << "Iter : " << solver.get_perfs().nb_iterations << std::endl
// << "Residuals : " << solver.get_perfs().residuals << std::endl;
if (cnt % 100 == 0)
{
scalar_t time = timestepper.get_total()/3600.0/24.0;
scalar_t sqrt_time = std::sqrt(time);
output_c_ca << time;
output_s_ca << time;
output_s_si << time;
output_loss_ca << time << "\t" << sqrt_time;
output_poro << time;
scalar_t amount_ca = 0.0;
for (index_t node: the_mesh->range_nodes())
{
amount_ca += variables->solid_concentration(node, id_ca, variables->displacement())
* the_mesh->get_volume_cell(node);
output_c_ca << "\t" << variables->aqueous_concentration(node, id_ca, variables->displacement());
output_s_ca << "\t" << variables->solid_concentration( node, id_ca, variables->displacement());
output_s_si << "\t" << variables->solid_concentration( node, id_si, variables->displacement());
output_poro << "\t" << variables->porosity(node);
}
output_poro << std::endl;
output_loss_ca << "\t" << (initial_ca - amount_ca)/1.75929e-2 << std::endl;
output_c_ca << std::endl;
output_s_ca << std::endl;
output_s_si << std::endl;
}
++cnt;
}
io::print_reactmicp_timer(&std::cout, solver.get_timer());
}
c4science · Help