Homec4science
Diffusion SpecMiCP / ReactMiCP (oxydoreduc) Browse

specmicp/examples/reactmicp/leaching.cpp
5246c20b3234oxydoreduc

Last Change
Raw File

leaching.cpp
View Options

#include <iostream>
#include "database/database.hpp"
#include "specmicp/reaction_path.hpp"
#include "reactmicp/mesh.hpp"
#include "reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp"
#include "reactmicp/systems/saturated_diffusion/boundary_conditions.hpp"
#include "reactmicp/systems/saturated_diffusion/reactive_transport_neutrality_solver.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_neutrality_parameters.hpp"
#include <fstream>
#include <Eigen/LU>
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::systems::siasaturated;
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);
}
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[-]"},
{"Al[3+]","Al(OH)4[-]"}
});
database.swap_components(swapping);
//std::vector<std::string> to_keep = {"H2O", "HO[-]", "Ca[2+]", "SiO(OH)3[-]"};
//database.keep_only_components(to_keep);
return database.get_database();
}
EquilibriumState get_initial_state(
std::shared_ptr<specmicp::database::DataContainer> database,
Vector& compo,
double mult)
{
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
compo = mult*compo;
const double m_c3s = compo(0);
const double m_c2s = compo(1);
const double m_c3a = compo(2);
const double m_gypsum = compo(3);
const double wc = 1.0;
const double m_water = wc*( (3*M_CaO+M_SiO2)*m_c3s
+ (2*M_CaO+M_SiO2)*m_c2s
+ (3*M_CaO+M_Al2O3)*m_c3a
+ (M_CaO + M_SO3)*m_gypsum
)*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)},
};
model->amount_minerals = {
{"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)},
{"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)},
{"C3A", specmicp::reaction_amount_t(m_c3a, 0.0)},
{"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0.0)},
};
model->minerals_to_keep = {
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2,am",
"Al(OH)3,am",
"Monosulfoaluminate",
"Straetlingite",
"Gypsum",
"Ettringite"
};
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 = 4;
driver.get_options().solver_options.factor_gradient_search_direction = 200;
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)
{
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
const double m_c3s = 1e-6;
const double m_c2s = 1e-6;
const double m_c3a = 1e-6;
const double m_gypsum = 1e-6;
const double m_water = 1;
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)},
{"SiO(OH)3[-]", specmicp::reaction_amount_t(0.0, 0.0)}
};
model->amount_minerals = {
{"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)},
{"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)},
{"C3A", specmicp::reaction_amount_t(m_c3a, 0.0)},
{"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0.0)},
//{"SiO2,am", specmicp::reaction_amount_t(0.0, 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();
}
double get_dt(
std::shared_ptr<SaturatedDiffusionTransportParameters> parameters,
std::shared_ptr<mesh::Mesh1D> themesh
)
{
const scalar_t max_d = parameters->diffusion_coefficient(1);
scalar_t dt = std::pow(themesh->get_dx(0), 2)/max_d/3;
if (dt > 0.01) return 0.01;
//return 10.0;
//return dt;
//if (dt<3600*24) dt=3600*24;
//if (dt>2*3600*24) dt=2*3600*24;
//return 900;
}
int leaching()
{
Vector compo_oxyde(4), compo_species(4);
compo_oxyde << 62.9, 20.6, 5.8, 3.1;
compo_from_oxyde(compo_oxyde, compo_species);
std::shared_ptr<database::DataContainer> thedatabase = set_database();
scalar_t radius = 3.5; //cm
scalar_t height = 8; //cm
index_t nb_nodes = 100;
mesh::Mesh1DPtr themesh = mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, height);
std::shared_ptr<SaturatedDiffusionTransportParameters> parameters =
std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e4*std::exp(9.95*0.15-29.08), 0.15);
parameters->diffusion_coefficient(0) = 2.2e-5;
parameters->porosity(0) = 1.0;
EquilibriumState initial_state = get_initial_state(thedatabase, compo_species, 0.5);
SIABoundaryConditions bcs(nb_nodes);
bcs.list_initial_states.push_back(initial_state);
bcs.list_initial_states.push_back(get_blank_state(thedatabase));
bcs.bs_types[0] = BoundaryConditionType::FixedComposition;
bcs.initial_states[0] = 1;
database::Database dbhandler(thedatabase);
index_t id_ca = dbhandler.component_label_to_id("Ca[2+]");
index_t id_si = dbhandler.component_label_to_id("SiO(OH)3[-]");
index_t id_oh = dbhandler.component_label_to_id("HO[-]");
std::ofstream output_ca, output_aqca, output_si, output_poro, output_ph;
output_ca.open("out_ca.dat");
output_aqca.open("out_aqca.dat");
output_si.open("out_si.dat");
output_ph.open("out_ph.dat");
output_poro.open("out_poro.ca");
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.get_variables().scale_volume();
solver.use_snia();
solver.use_sia(100);
solver.get_options().residual_tolerance = 1e-3;
solver.get_options().threshold_update_disable_speciation = 1e-16;
solver.get_options().threshold_update_transport = 1e-10;
// Transport
solver.get_options().transport_options.residuals_tolerance = 1e-6;
solver.get_options().transport_options.quasi_newton = 5;
//solver.get_options().transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
// Speciation
solver.get_options().speciation_options.solver_options.use_crashing = false;
solver.get_options().speciation_options.solver_options.use_scaling = false;
solver.get_options().speciation_options.solver_options.max_iter = 100;
solver.get_options().speciation_options.solver_options.maxstep = 50;
solver.get_options().speciation_options.solver_options.non_monotone_linesearch = true;
std::cout << "id SiO2,am : " << dbhandler.mineral_label_to_id("SiO2,am") << std::endl;
//thedatabase->stability_mineral(
// dbhandler.mineral_label_to_id("SiO2,am")) = database::MineralStabilityClass::cannot_dissolve;
//double dt = 5e4;
double total=0 ;
uindex_t k = 1;
while (total < 100*24*3600)
{
double dt = get_dt(parameters, themesh);
total += dt;
std::cout << " iterations : " << k << " - dt : " << dt << 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.total_mobile_concentrations() << std::endl;
// compute porosity
for (index_t node: themesh->range_nodes())
{
double volume = variables.volume_minerals(node);
double cell_volume = themesh->get_volume_cell(node);
specmicp_assert(volume < cell_volume);
double porosity = (cell_volume - volume*1e6) / cell_volume;
parameters->porosity(node) = porosity;
parameters->diffusion_coefficient(node) =
porosity>0.92?1e4*std::exp(9.95*0.92-29.08):1e4*std::exp(9.95*porosity-29.08);
//std::cout << parameters->diffusion_coefficient(node);
}
parameters->diffusion_coefficient(0) =1e4*std::exp(9.95*0.92-29.08); // parameters->diffusion_coefficient(1);
//std::cout << variables.speciation_variables() << std::endl;
if (k==1 or (k < 1000 and k% 100 == 0) or k % 1000 == 0)
{
output_ca << total << " " << total/(3600.0*24.0);
output_aqca << total << " " << total/(3600.0*24.0);
output_si << total << " " << total/(3600.0*24.0);
output_poro << total << " " << total/(3600.0*24.0);
output_ph << total << " " << total/(3600.0*24.0);
for(index_t node: themesh->range_nodes())
{
output_ca << " " << variables.nodal_component_update_total_amount(node, id_ca);
output_aqca << " " << variables.total_mobile_concentration(node, id_ca);
output_si << " " << variables.nodal_component_update_total_amount(node, id_si);
output_ph << " " <<
14 + variables.log_gamma_component(node, id_oh)
+ variables.log_component_concentration(node, id_oh);
double volume = variables.volume_minerals(node);
double cell_volume = themesh->get_volume_cell(node);
output_poro << " " << (cell_volume - volume*1e6) / cell_volume;
}
output_ca << std::endl;
output_aqca << std::endl;
output_si << std::endl;
output_poro << std::endl;
output_ph << std::endl;
for (index_t mineral: thedatabase->range_mineral())
{
mineral_out[mineral] << total << " " << total/(3600.0*24.0);
for(index_t node: themesh->range_nodes())
{
mineral_out[mineral] << " " << variables.mineral_amount(node, mineral);
}
mineral_out[mineral] << std::endl;
}
}
++k;
}
output_ca.close();
output_aqca.close();
output_si.close();
output_poro.close();
for (index_t mineral: thedatabase->range_mineral())
{
mineral_out[mineral].close();
}
return EXIT_SUCCESS;
}
int leaching_neutrality()
{
std::shared_ptr<database::DataContainer> thedatabase = set_database();
Vector compo_oxyde(4), compo_species(4);
compo_oxyde << 62.9, 20.6, 5.8, 3.1;
compo_from_oxyde(compo_oxyde, compo_species);
scalar_t radius = 3.5; //cm
scalar_t height = 8; //cm
index_t nb_nodes = 36;
database::Database dbhandler(thedatabase);
//mesh::Mesh1DPtr themesh = mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, height);
mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 5);
//auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-8, 0.2);
std::shared_ptr<SaturatedNeutralityDiffusionTransportParameters> parameters =
std::make_shared<SaturatedNeutralityDiffusionTransportParameters>(
nb_nodes, // nb_nodes,
thedatabase->nb_component, // nb_component,
thedatabase->nb_aqueous, // nb_aqueous,
2e-5, // the_diffusion_coefficient,
0.25, // porosity
0.5e-3 //reduction_factor
);
//parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("H[+]")) = 9e-5;
//parameters->diff_coeff_component(dbhandler.component_label_to_id("HO[-]")) = 9e-5;
parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("Si(OH)4")) = 1e-7;
parameters->diff_coeff_component(dbhandler.component_label_to_id("SiO(OH)3[-]")) = 1e-7;
parameters->porosity(0) = 1.0;
parameters->reduction_factor(0) = 1.0;
EquilibriumState initial_state = get_initial_state(thedatabase, compo_species, 0.5);
SIABoundaryConditions bcs(nb_nodes);
bcs.list_initial_states.push_back(initial_state);
bcs.list_initial_states.push_back(get_blank_state(thedatabase));
bcs.bs_types[0] = BoundaryConditionType::FixedComposition;
bcs.initial_states[0] = 1;
index_t id_ca = dbhandler.component_label_to_id("Ca[2+]");
index_t id_si = dbhandler.component_label_to_id("SiO(OH)3[-]");
std::ofstream output_ca, output_aqca, output_si, output_poro;
output_ca.open("out_ca.dat");
output_aqca.open("out_aqca.dat");
output_si.open("out_si.dat");
output_poro.open("out_poro.ca");
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");
}
SIASaturatedReactiveTransportNeutralitySolver solver(themesh, thedatabase, parameters);
solver.apply_boundary_conditions(bcs);
solver.get_variables().scale_volume();
solver.use_sia(100);
solver.get_options().residual_tolerance = 1e-3;
solver.get_options().threshold_update_disable_speciation = 1e-16;
solver.get_options().threshold_update_transport = 1e-6;
// Transport
solver.get_options().transport_options.residuals_tolerance = 1e-5;
solver.get_options().transport_options.maximum_iterations = 100;
solver.get_options().transport_options.quasi_newton = 2;
solver.get_options().transport_options.maximum_step_length = 1000;
solver.get_options().transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
solver.get_options().transport_options.max_iterations_at_max_length = 100;
// Speciation
solver.get_options().speciation_options.solver_options.use_crashing = false;
solver.get_options().speciation_options.solver_options.use_scaling = false;
solver.get_options().speciation_options.solver_options.max_iter = 100;
solver.get_options().speciation_options.solver_options.maxstep = 50;
std::cout << "id SiO2,am : " << dbhandler.mineral_label_to_id("SiO2,am") << std::endl;
//thedatabase->stability_mineral(
// dbhandler.mineral_label_to_id("SiO2,am")) = database::MineralStabilityClass::cannot_dissolve;
//double dt = 5e4;
double total=0 ;
uindex_t k = 1;
std::cout << "Basis :";
for (index_t component: thedatabase->range_component())
{
std:: cout << " " << thedatabase->labels_basis[component];
}
std::cout << std::endl;
std::cout << "Solid :";
for (index_t mineral: thedatabase->range_mineral())
{
std:: cout << " " << thedatabase->labels_minerals[mineral];
}
std::cout << std::endl;
while (total < 10*24*3600)
{
double dt = 3600/2;
total += dt;
std::cout << " iterations : " << k << " - dt : " << dt << 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.total_mobile_concentrations() << std::endl;
// compute porosity
for (index_t node: themesh->range_nodes())
{
double volume = variables.volume_minerals(node);
double cell_volume = themesh->get_volume_cell(node);
specmicp_assert(volume < cell_volume);
double porosity = (cell_volume - volume*1e6) / cell_volume;
parameters->porosity(node) = porosity;
parameters->reduction_factor(node) =
porosity>0.6?1e4*std::exp(9.95*0.6-29.08)/2e-5:1e4*std::exp(9.95*porosity-29.08)/2e-5;
//std::cout << parameters->diffusion_coefficient(node);
}
parameters->reduction_factor(0) = parameters->reduction_factor(1);
//std::cout << variables.speciation_variables() << std::endl;
if (k < 2 or k % 100 == 0)
{
output_ca << total << " " << total/(3600.0*24.0);
output_aqca << total << " " << total/(3600.0*24.0);
output_si << total << " " << total/(3600.0*24.0);
output_poro << total << " " << total/(3600.0*24.0);
for(index_t node: themesh->range_nodes())
{
output_ca << " " << variables.nodal_component_update_total_amount(node, id_ca);
output_aqca << " " << variables.total_mobile_concentration(node, id_ca);
output_si << " " << variables.nodal_component_update_total_amount(node, id_si);
double volume = variables.volume_minerals(node);
double cell_volume = themesh->get_volume_cell(node);
output_poro << " " << (cell_volume - volume*1e6) / cell_volume;
}
output_ca << std::endl;
output_aqca << std::endl;
output_si << std::endl;
output_poro << std::endl;
for (index_t mineral: thedatabase->range_mineral())
{
mineral_out[mineral] << total << " " << total/(3600.0*24.0);
for(index_t node: themesh->range_nodes())
{
mineral_out[mineral] << " " << variables.mineral_amount(node, mineral);
}
mineral_out[mineral] << std::endl;
}
}
++k;
}
output_ca.close();
output_aqca.close();
output_si.close();
output_poro.close();
for (index_t mineral: thedatabase->range_mineral())
{
mineral_out[mineral].close();
}
return EXIT_SUCCESS;
}
int main()
{
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
specmicp::logger::ErrFile::stream() = &std::cerr;
return leaching();
//return leaching_neutrality();
}
c4science · Help