Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81416503
leaching.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
Thu, Sep 5, 23:17
Size
19 KB
Mime Type
text/x-c
Expires
Sat, Sep 7, 23:17 (2 d)
Engine
blob
Format
Raw Data
Handle
20563547
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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-7;
const double m_c2s = 1e-7;
const double m_c3a = 1e-7;
const double m_gypsum = 1e-7;
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.0001, 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;
//if (dt<3600*24) dt=3600*24;
//if (dt>2*3600*24) dt=2*3600*24;
return 3600*24/5;
}
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 = 36;
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.25-29.08), 0.25);
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[-]");
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");
}
SIASaturatedReactiveTransportSolver 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-10;
// Transport
solver.get_options().transport_options.residuals_tolerance = 1e-5;
solver.get_options().transport_options.quasi_newton = 3;
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;
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 < 365*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.6?1e4*std::exp(9.95*0.6-29.08):1e4*std::exp(9.95*porosity-29.08);
//std::cout << parameters->diffusion_coefficient(node);
}
parameters->diffusion_coefficient(0) = parameters->diffusion_coefficient(1);
//std::cout << variables.speciation_variables() << std::endl;
if (k < 10 or k % 5 == 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 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 < 10 or k % 5 == 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_neutrality();
}
Event Timeline
Log In to Comment