Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F63518419
carbonation.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
Mon, May 20, 17:01
Size
20 KB
Mime Type
text/x-c
Expires
Wed, May 22, 17:01 (2 d)
Engine
blob
Format
Raw Data
Handle
17745593
Attached To
rSPECMICP SpecMiCP / ReactMiCP
carbonation.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 "utils/moving_average.hpp"
#include <fstream>
#include <Eigen/LU>
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(-1e-5, 0.0)},
{"HCO3[-]", specmicp::reaction_amount_t(1e-5, 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;
driver.get_options().system_options.charge_keeper = 1;
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(const EquilibriumState& healthy_cp)
{
specmicp::Vector totconc;
specmicp::RawDatabasePtr thedatabase = healthy_cp.get_database();
specmicp::Database dbhandler(thedatabase);
// healthy_cp.total_concentrations(totconc);
// index_t id_hco3 = dbhandler.component_label_to_id("HCO3[-]");
// index_t id_ca = dbhandler.component_label_to_id("Ca[2+]");
index_t id_oh = dbhandler.component_label_to_id("HO[-]");
// index_t id_ch = dbhandler.mineral_label_to_id("Portlandite");
// index_t id_csh = dbhandler.mineral_label_to_id("CSH,jennite");
// scalar_t addition =
// 2.5*(
// thedatabase->nu_mineral(id_ch, id_ca)*healthy_cp.moles_mineral(id_ch)
// + thedatabase->nu_mineral(id_csh, id_ca)*healthy_cp.moles_mineral(id_csh)
// );
// totconc(id_hco3) += addition;
// totconc(id_oh) -= addition;
// specmicp::Vector variables(healthy_cp.main_variables());
// specmicp::ReducedSystemSolver solver(thedatabase, totconc, healthy_cp);
// solver.get_options().system_options.charge_keeper = id_oh;
// specmicp::micpsolver::MiCPPerformance perf = solver.solve(variables);
// if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized)
// throw std::runtime_error("Error - blank system cannot be solved");
// return solver.get_solution(variables);
std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
//const double m_c3s = 1e-8;
const double m_c2s = 1e-6;
const double m_water = 1.0;
model->amount_aqueous = {
{"H2O", specmicp::reaction_amount_t(m_water/thedatabase->molar_mass_basis_si(0), 0.0)},
// {"HO[-]", specmicp::reaction_amount_t(-0.2, 0.0)},
// {"SiO(OH)3[-]", specmicp::reaction_amount_t(0.0, 0.0)},
// {"HCO3[-]", specmicp::reaction_amount_t(0.2, 0.0)}
};
model->amount_minerals = {
//{"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)},
{"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)},
{"SiO2,am", specmicp::reaction_amount_t(1e-2, 0.0)},
{"Calcite", specmicp::reaction_amount_t(1e-2, 0.0)}
};
Eigen::VectorXd x(thedatabase->nb_component+thedatabase->nb_mineral);
x.setZero();
x(0) = 1.0;
x.segment(1, thedatabase->nb_component-1).setConstant(-6);
specmicp::ReactionPathDriver driver(model, thedatabase, 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().system_options.charge_keeper = id_oh;
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;
return solution;
}
int carbonation()
{
//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 = 100;
utils::ExponentialMovingAverage averager1(0.001, 5);
utils::ExponentialMovingAverage averager2(0.01, 5);
utils::ExponentialMovingAverage averager3(0.1, 5);
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.25-29.08),
0.25);
parameters->diffusion_coefficient(0) = 2.2e-5; // cm^2/s
parameters->porosity(0) = 1.0;
EquilibriumState initial_state = get_initial_state(thedatabase, 0.5);
SIABoundaryConditions bcs(nb_nodes);
bcs.list_initial_states.push_back(initial_state);
bcs.list_initial_states.push_back(get_blank_state(initial_state));
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[-]");
index_t id_co2 = dbhandler.component_label_to_id("HCO3[-]");
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;
std::ofstream output_co2, output_aqco2, output_sa_co2;
std::ofstream output_iter;
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");
output_iter.open("nb_iter.dat");
output_co2.open("out_co2.dat");
output_aqco2.open("out_aqco2.dat");
output_sa_co2.open("out_sa_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.get_variables().scale_volume();
//solver.use_snia();
solver.use_sia(500);
solver.get_options().residual_tolerance = 1e-2;
solver.get_options().absolute_residual_tolerance = 1e-8;
solver.get_options().threshold_update_disable_speciation = 1e-20;
solver.get_options().threshold_update_transport = 1e-10;
solver.get_options().good_enough_transport_residual_tolerance = 1.0;
solver.get_options().relative_update = 1e-6;
solver.get_options().use_consistent_velocity = true;
solver.get_options().use_previous_flux = false;
// Transport
solver.get_options().transport_options.residuals_tolerance = 1e-8;
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 = 3;
solver.get_options().transport_options.sparse_solver = SparseSolver::SparseLU;
solver.get_options().transport_options.maximum_step_length = 1e4;
//solver.get_options().transport_options.alpha = 0.0;
solver.get_options().transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
// Speciation
solver.get_options().speciation_options.system_options.charge_keeper = 1;
solver.get_options().speciation_options.system_options.conservation_water = true;
solver.get_options().speciation_options.system_options.non_ideality = true;
solver.get_options().speciation_options.solver_options.fvectol = 1e-12;
solver.get_options().speciation_options.system_options.non_ideality_tolerance = 1e-14;
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 = 40;
solver.get_options().speciation_options.solver_options.non_monotone_linesearch = true;
solver.get_options().speciation_options.solver_options.threshold_cycling_linesearch = 1e-8;
solver.get_options().speciation_options.solver_options.steptol = 1e-16;
solver.get_options().speciation_options.solver_options.threshold_stationary_point = 1e-4;
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 = 10.0;
scalar_t dt_min = 1.0;
scalar_t dt_max = 10.0;
uindex_t max_iter;
while (total < 365*24*3600)
{
total += dt;
std::cout << " iterations : " << k << " - dt : " << dt << " - total : " << total/24/3600 << std::endl;
//start_again:
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)
{
if (false and retcode == SIASaturatedReactiveTransportSolverReturnCode::MaxIterations)
{
max_iter += 1;
ERROR << "Maximum iterations";
}
else
{
std::cout << "Current residual : " << solver.get_perfs().current_residual << std::endl;
SIASaturatedVariables& variables = solver.get_variables();
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;
output_ca.close();
output_aqca.close();
output_si.close();
output_poro.close();
output_saca.close();
output_sasi.close();
throw std::runtime_error("Failed to converge : "+std::to_string((int) retcode));
}
}
//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) =
// std::max(porosity>0.92?2.219e-5:1e4*std::exp(9.95*porosity-29.08), 2e-8);
}
std::cout << variables.component_concentration(1, id_ca) << std::endl;
parameters->diffusion_coefficient(0) = 2.219e-5; // parameters->diffusion_coefficient(1);
//if (solver.get_perfs().nb_iterations > 15) dt = 0.95*dt;
//else if (solver.get_perfs().nb_iterations < 2) dt = 1.05*dt;
output_iter << k
<< " " << solver.get_perfs().nb_iterations
<< " " << solver.get_perfs().current_residual
<< " " << averager1.add_point(solver.get_perfs().nb_iterations)
<< " " << averager2.add_point(solver.get_perfs().nb_iterations)
<< " " << averager3.add_point(solver.get_perfs().nb_iterations)
<< std::endl;
//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_saca << total << " " << total/(3600.0*24.0);
output_si << total << " " << total/(3600.0*24.0);
output_sasi << total << " " << total/(3600.0*24.0);
output_co2 << total << " " << total/(3600.0*24.0);
output_aqco2 << total << " " << total/(3600.0*24.0);
output_sa_co2 << 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_saca << " " << variables.component_concentration_in_minerals(node, id_ca);
output_sasi << " " << variables.component_concentration_in_minerals(node, id_si);
output_si << " " << variables.nodal_component_update_total_amount(node, id_si);
output_co2 << " " << variables.nodal_component_update_total_amount(node, id_co2);
output_aqco2 << " " << variables.total_mobile_concentration(node, id_co2);
output_sa_co2 << " " << variables.component_concentration_in_minerals(node, id_co2);
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_saca << std::endl;
output_si << std::endl;
output_sasi << std::endl;
output_co2 << std::endl;
output_aqco2 << std::endl;
output_sa_co2 << 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;
if ((averager3.current_value() - averager1.current_value()) > 10)
{
dt = std::max(dt*0.95,dt_min);
averager1.reset(averager3.current_value());
}
if ((averager2.current_value() - averager3.current_value()) > 10)
{
dt = std::min(dt*1.5, dt_max);
averager2.reset(averager3.current_value());
}
}
output_iter.close();
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();
}
std::cout << "Max iterations : " << max_iter << std::endl;
return EXIT_SUCCESS;
}
int main()
{
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
specmicp::logger::ErrFile::stream() = &std::cerr;
carbonation();
}
Event Timeline
Log In to Comment