Homec4science
Diffusion SpecMiCP / ReactMiCP (fixpoint) Browse

specmicp/examples/reactmicp/carbonation.cpp
edda4626caa0fixpoint

Last Change
Raw File

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();
}
c4science · Help