Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/examples/reactmicp/carbonation.cpp
6d08d8744732reactivetransport

Last Change
Raw File

carbonation.cpp
View Options

#include <fstream>
#include "database/database.hpp"
#include "specmicp/reaction_path.hpp"
#include "reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp"
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::systems::siasaturated;
std::shared_ptr<database::DataContainer> set_database() {
specmicp::database::Database database("../data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
});
database.swap_components(swapping);
//std::vector<std::string> to_keep = {"H2O", "HO[-]", "Ca[2+]", "SiO(OH)3[-]", "HCO3[-]", "Al(OH)4[-]", "SO2"};
//database.keep_only_components(to_keep);
return database.get_database();
}
EquilibriumState get_initial_state(std::shared_ptr<specmicp::database::DataContainer> database, double mult)
{
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
const double m_c3s = mult*0.7;
const double m_c2s = mult*0.3;
const double wc = 1.0;
const double m_water = wc*((3*56.08+60.08)*m_c3s + (2*56.08+60.08)*m_c2s)*1e-3;
model->amount_aqueous = {
{"H2O", specmicp::reaction_amount_t(m_water/database->molar_mass_basis_si(0), 0.0)},
{"HO[-]", specmicp::reaction_amount_t(0.0, 0.0)},
{"HCO3[-]", specmicp::reaction_amount_t(1e-9, 0.0)}
};
model->amount_minerals = {
{"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)},
{"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)}
};
Eigen::VectorXd x;
specmicp::ReactionPathDriver driver(model, database, x);
//driver.get_options().solver_options.penalization_factor =1.0;
driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
//driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
driver.get_options().solver_options.use_scaling = false;
driver.get_options().solver_options.max_factorization_step = 1;
driver.get_options().solver_options.factor_gradient_search_direction = 50;
driver.get_options().solver_options.maxstep = 50;
driver.get_options().solver_options.maxiter_maxstep = 50;
driver.get_options().solver_options.max_iter = 100;
driver.get_options().allow_restart = true;
specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x);
if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized)
throw std::runtime_error("Error - system cannot be solved");
return driver.get_current_solution();
}
EquilibriumState get_blank_state(std::shared_ptr<specmicp::database::DataContainer> database, double mult)
{
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
const double m_c3s = mult*0.7;
const double m_c2s = mult*0.3;
const double wc = 1.0;
const double m_water = wc*((3*56.08+60.08)*m_c3s + (2*56.08+60.08)*m_c2s)*1e-3;
const double t_hco3 = mult*2.6;
model->amount_aqueous = {
{"H2O", specmicp::reaction_amount_t(m_water/database->molar_mass_basis_si(0), 0.0)},
{"HO[-]", specmicp::reaction_amount_t(-t_hco3, 0.0)},
{"HCO3[-]", specmicp::reaction_amount_t(t_hco3, 0.0)}
};
model->amount_minerals = {
{"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)},
{"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)}
};
Eigen::VectorXd x;
specmicp::ReactionPathDriver driver(model, database, x);
//driver.get_options().solver_options.penalization_factor =1.0;
driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
//driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
driver.get_options().solver_options.use_scaling = false;
driver.get_options().solver_options.max_factorization_step = 1;
driver.get_options().solver_options.factor_gradient_search_direction = 50;
driver.get_options().solver_options.maxstep = 50;
driver.get_options().solver_options.maxiter_maxstep = 50;
driver.get_options().solver_options.max_iter = 100;
driver.get_options().allow_restart = true;
specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x);
if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized)
throw std::runtime_error("Error - blank system cannot be solved");
return driver.get_current_solution();
}
int main()
{
std::shared_ptr<database::DataContainer> thedatabase = set_database();
database::Database dbhandler(thedatabase);
int id_ca = dbhandler.component_label_to_id("Ca[2+]");
int id_co2 = dbhandler.component_label_to_id("HCO3[-]");
int nb_nodes = 50;
std::shared_ptr<mesh::Mesh1D> themesh = std::make_shared<mesh::Mesh1D>(nb_nodes-1, 0.1, 2.5);
auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(1e-4, 0.2);
EquilibriumState initial_state = get_initial_state(thedatabase, 0.02);
// scaling {
/*
double volume = initial_state.volume_minerals();
std::cout << volume << std::endl;
volume /= (1 - parameters->porosity(1) );
double scale = themesh->cell_volume(1)*1e-6 / volume;
initial_state.scale_condensed(scale);
std::cout << initial_state.volume_minerals() << std::endl;
*/
// }
SIABoundaryConditions bcs(nb_nodes);
bcs.list_initial_states.push_back(initial_state);
bcs.list_initial_states.push_back(get_blank_state(thedatabase, 0.02));
bcs.bs_types[0] = BoundaryConditionType::FixedComposition;
bcs.initial_states[0] = 1;
std::ofstream output_ca;
output_ca.open("out_ca.dat");
std::ofstream output_co2;
output_co2.open("out_co2.dat");
std::vector<std::ofstream> mineral_out(thedatabase->nb_mineral);
for (int mineral: thedatabase->range_mineral())
{
mineral_out[mineral].open("out_mineral_"+std::to_string(mineral)+".dat");
}
SIASaturatedReactiveTransportSolver solver(themesh, thedatabase, parameters);
solver.apply_boundary_conditions(bcs);
solver.use_sia(100);
double dt = 100;
double total=0 ;
for (int k=0; k<10000; ++k)
{
std::cout << " iterations : " << k << std::endl;
SIASaturatedReactiveTransportSolverReturnCode retcode = solver.solve_timestep(dt);
if (retcode != SIASaturatedReactiveTransportSolverReturnCode::Success) throw std::runtime_error("Failed to converge : "+std::to_string((int) retcode));
SIASaturatedVariables& variables = solver.get_variables();
//std::cout << variables.speciation_variables() << std::endl;
if (k % 5 == 0)
{
output_ca << total << " " << total/(3600.0*24.0);
output_co2 << total << " " << total/(3600.0*24.0);
for(int node: themesh->range_nodes())
{
output_ca << " " << variables.nodal_component_update_total_amount(node, id_ca);
output_co2 << " " << variables.nodal_component_update_total_amount(node, id_co2);
}
output_ca << std::endl;
output_co2 << std::endl;
for (int mineral: thedatabase->range_mineral())
{
mineral_out[mineral] << total << " " << total/(3600.0*24.0);
for(int node: themesh->range_nodes())
{
mineral_out[mineral] << " " << variables.mineral_amount(node, mineral);
}
mineral_out[mineral] << std::endl;
}
}
total+=dt;
}
output_ca.close();
output_co2.close();
for (int mineral: thedatabase->range_mineral())
{
mineral_out[mineral].close();
}
}
c4science · Help