Homec4science
Diffusion SpecMiCP / ReactMiCP (bin/unsaturated) Browse

specmicp/examples/reactmicp/carbonation.cpp
5246c20b3234bin/unsaturated

Last Change
Raw File

carbonation.cpp
View Options

#include <fstream>
#include "database/database.hpp"
#include "specmicp/reaction_path.hpp"
#include "reactmicp/meshes/uniform_mesh1d.hpp"
#include "reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp"
#include "reactmicp/systems/saturated_diffusion/boundary_conditions.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_program.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp"
#include <iostream>
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(0.0, 0.0)},
{"HCO3[-]", specmicp::reaction_amount_t(1e-11, 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().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 =
3.0*(
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().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);
}
int main()
{
std::shared_ptr<database::DataContainer> thedatabase = set_database();
database::Database dbhandler(thedatabase);
int nb_nodes = 100;
std::shared_ptr<mesh::Mesh1D> themesh = std::make_shared<mesh::UniformMesh1D>(nb_nodes, 0.1, 5.0);
auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-5, 0.2);
specmicp::Vector totconc;
EquilibriumState initial_state = get_initial_state(thedatabase, 0.02);
int id_ca = dbhandler.component_label_to_id("Ca[2+]");
int id_co2 = dbhandler.component_label_to_id("HCO3[-]");
std::cout << id_co2 << std::endl;
std::cout << thedatabase->nu_aqueous.col(id_co2) << std::endl;
initial_state.total_concentrations(totconc);
std::cout << "tot co2 : " << totconc(id_co2) << std::endl;
for (index_t aqueous=0; aqueous<thedatabase->nb_aqueous; ++aqueous)
{
std::cout << thedatabase->labels_aqueous[aqueous] << " - "
<< initial_state.molality_secondary(aqueous) << std::endl;
}
// scaling {
double volume = initial_state.volume_minerals();
std::cout << volume << std::endl;
volume /= (1 - parameters->porosity(1) );
double scale = themesh->get_volume_cell(1)*1e-6 / volume;
initial_state.scale_condensed(scale);
std::cout << initial_state.volume_minerals() << std::endl;
// }
initial_state.total_concentrations(totconc);
std::cout << "tot co2 : " << totconc(id_co2) << std::endl;
EquilibriumState bc_state = get_blank_state(initial_state);
SIABoundaryConditions bcs(nb_nodes);
bcs.list_initial_states.push_back(initial_state);
bcs.list_initial_states.push_back(bc_state);
bcs.bs_types[0] = BoundaryConditionType::FixedComposition;
bcs.initial_states[0] = 1;
std::ofstream output_ca;
output_ca.open("out_ca.dat");
std::ofstream output_co2;
output_co2.open("out_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_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 = 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;
solver.use_sia(30);
double dt = 50;
double total=0 ;
int k=0;
while(total < 3600*24*100)
{
std::cout << " iterations : " << k << 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.speciation_variables() << std::endl;
if (k==1 or k % 100 == 0)
{
output_ca << total << " " << total/(3600.0*24.0);
output_co2 << total << " " << total/(3600.0*24.0);
for(int node: themesh->range_nodes())
{
output_ca << " " << variables.nodal_component_update_total_amount(node, id_ca);
output_co2 << " " << variables.nodal_component_update_total_amount(node, id_co2);
}
output_ca << std::endl;
output_co2 << std::endl;
for (int mineral: thedatabase->range_mineral())
{
mineral_out[mineral] << total << " " << total/(3600.0*24.0);
for(int node: themesh->range_nodes())
{
mineral_out[mineral] << " " << variables.mineral_amount(node, mineral);
}
mineral_out[mineral] << std::endl;
}
}
total+=dt;
++k;
}
output_ca.close();
output_co2.close();
for (int mineral: thedatabase->range_mineral())
{
mineral_out[mineral].close();
}
}
c4science · Help