Homec4science
Diffusion SpecMiCP / ReactMiCP (configuration) Browse

specmicp/examples/reactmicp/leaching.cpp
bad60df63f6fconfiguration

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-8;
//const double m_c2s = 1e-8;
//const double m_c3a = 10e-8;
//const double m_gypsum = 10e-8;
//const double m_water = 1.0;
//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.0001, 0.0)}
//};
//Eigen::VectorXd x(database->nb_component+database->nb_mineral);
//x.setZero();
//x(0) = 1.0;
//x.segment(1, database->nb_component-1).setConstant(-6);
//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");
//
//EquilibriumState solution = driver.get_current_solution();
//std::cout << solution.main_variables() << std::endl;
//solution.remove_minerals();
Eigen::VectorXd variables(database->nb_component+database->nb_mineral);
variables.setZero();
variables(0) = 1.0;
for (index_t component = 1; component < database->nb_component; ++component)
variables(component) = -11;
scalar_t ionic_strength = 0.0;
Eigen::VectorXd secondary(database->nb_aqueous);
Eigen::VectorXd loggamma(database->nb_component+database->nb_aqueous);
secondary.setZero();
loggamma.setZero();
EquilibriumState solution(variables, secondary, loggamma, ionic_strength, database);
return 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 > 5.0) return 5.0;
//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.0; //cm
scalar_t dx = 0.005;
index_t additional_nodes = 1;
radius = radius + additional_nodes*dx;
index_t nb_nodes = 50;
mesh::Mesh1DPtr themesh = mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, dx, 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; // cm^2/s
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;
for (index_t node =0; node<additional_nodes; ++node)
bcs.initial_states[node] = 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[-]");
for (index_t mineral: thedatabase->range_mineral())
{
std::cout << thedatabase->labels_minerals[mineral] << std::endl;
}
std::cout << "Volume first cell " << themesh->get_volume_cell(additional_nodes) << std::endl;
//return EXIT_SUCCESS;
std::ofstream output_ca, output_aqca, output_si, output_poro, output_ph;
std::ofstream output_saca, output_sasi;
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.dat");
output_saca.open("out_sa_ca.dat");
output_sasi.open("out_sa_si.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.get_variables().scale_volume();
//solver.use_snia();
solver.use_sia(100);
solver.get_options().residual_tolerance = 1e-6;
solver.get_options().threshold_update_disable_speciation = 1e-16;
solver.get_options().threshold_update_transport = 1e-9;
solver.get_options().good_enough_transport_residual_tolerance = 1e-3;
solver.get_options().relative_update = 1e-4;
// Transport
solver.get_options().transport_options.residuals_tolerance = 1e-12;
solver.get_options().transport_options.maximum_iterations = 500;
solver.get_options().transport_options.step_tolerance = 1e-8;
solver.get_options().transport_options.threshold_stationary_point = 1e-5;
solver.get_options().transport_options.quasi_newton = 5;
solver.get_options().transport_options.sparse_solver = SparseSolver::GMRES;
//solver.get_options().transport_options.alpha = 0.0;
solver.get_options().transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
// Speciation
solver.get_options().speciation_options.solver_options.fvectol = 1e-10;
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 = 50;
solver.get_options().speciation_options.solver_options.maxstep = 50;
solver.get_options().speciation_options.solver_options.non_monotone_linesearch = true;
solver.get_options().speciation_options.solver_options.threshold_cycling_linesearch = 1e-8;
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;
SIASaturatedVariables& first_variables = solver.get_variables();
std::cout << first_variables.total_mobile_concentrations() << std::endl;
//double dt = 5e4;
double total=0 ;
uindex_t k = 0;
scalar_t dt = 20.0;
while (total < 365*24*3600)
{
total += dt;
std::cout << " iterations : " << k << " - dt : " << dt << std::endl;
SIASaturatedReactiveTransportSolverReturnCode retcode = solver.solve_timestep(dt);
std::cout << " Return code : " << (int) retcode
<< " - Nb iter : " << solver.get_perfs().nb_iterations
<< std::endl
<< " - Nb iter transport : " << solver.get_perfs().nb_iterations_transport
<< std::endl
<< " - Nb iter chemistry : " << solver.get_perfs().nb_iterations_chemistry
<< std::endl;
if (retcode != SIASaturatedReactiveTransportSolverReturnCode::Success)
{
std::cout << "Current residual : " << solver.get_perfs().current_residual << std::endl;
throw std::runtime_error("Failed to converge : "+std::to_string((int) retcode));
//dt = dt/1.5;
//goto start_again;
}
//if (solver.get_perfs().nb_iterations <= 10) dt*=1.1;
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?2.219e-5:1e4*std::exp(9.95*porosity-29.08);
}
std::cout << variables.component_concentration(1, id_ca) << std::endl;
parameters->diffusion_coefficient(0) = 2.219e-5; // 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);
output_saca << total << " " << total/(3600.0*24.0);
output_sasi << 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_saca << " " << variables.component_concentration_in_minerals(node, id_ca);
output_sasi << " " << variables.component_concentration_in_minerals(node, id_si);
}
output_ca << std::endl;
output_aqca << std::endl;
output_si << std::endl;
output_poro << std::endl;
output_ph << std::endl;
output_saca << std::endl;
output_sasi << 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();
output_saca.close();
output_sasi.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