Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F120843722
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, Jul 7, 11:34
Size
19 KB
Mime Type
text/x-c++
Expires
Wed, Jul 9, 11:34 (2 d)
Engine
blob
Format
Raw Data
Handle
27258712
Attached To
rSPECMICP SpecMiCP / ReactMiCP
carbonation.cpp
View Options
#include <Eigen/Core>
#include <Eigen/QR>
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "reactmicp/systems/saturated_react/variables.hpp"
#include "reactmicp/systems/saturated_react/equilibrium_stagger.hpp"
#include "reactmicp/systems/saturated_react/transport_stagger.hpp"
#include "reactmicp/systems/saturated_react/init_variables.hpp"
#include "reactmicp/solver/reactive_transport_solver.hpp"
#include "reactmicp/solver/reactive_transport_solver_structs.hpp"
#include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp"
#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
#include "reactmicp/solver/timestepper.hpp"
#include "physics/laws.hpp"
#include "dfpm/mesh.hpp"
#include "utils/log.hpp"
#include "utils/timer.hpp"
#include "utils/io/reactive_transport.hpp"
#include <iostream>
#include <fstream>
using namespace specmicp;
using namespace reactmicp;
using namespace reactmicp::systems;
using namespace specmicp::reactmicp::systems::satdiff;
using namespace reactmicp::solver;
AdimensionalSystemSolution get_cement_initial_compo(
RawDatabasePtr the_database,
const AdimensionalSystemConstraints& constraints,
const AdimensionalSystemSolverOptions& options
);
AdimensionalSystemConstraints get_cement_initial_constraints(
RawDatabasePtr the_database,
const units::UnitsSet& units_set);
AdimensionalSystemSolution get_bc_initial_compo(
RawDatabasePtr the_database,
const AdimensionalSystemConstraints& constraints,
const AdimensionalSystemSolverOptions& options
);
AdimensionalSystemConstraints get_bc_initial_constraints(
RawDatabasePtr the_database,
const units::UnitsSet& units_set);
RawDatabasePtr get_database(const std::string& path_to_database);
Vector initial_compo(const Vector& publi);
AdimensionalSystemSolverOptions get_specmicp_options();
mesh::Mesh1DPtr get_mesh(scalar_t dx);
void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options);
void set_reactive_solver_options(reactmicp::solver::ReactiveTransportOptions& reactive_options);
void print_minerals_profile(
RawDatabasePtr the_database,
SaturatedVariablesPtr variables,
mesh::Mesh1DPtr the_mesh,
const std::string& filepath
);
class CarboUspcaling: public solver::UpscalingStaggerBase
{
public:
static constexpr scalar_t de_0 = 1e-11;
static constexpr scalar_t factor = 1e4;
static constexpr scalar_t porosity_r = 0.02;
CarboUspcaling(RawDatabasePtr the_database,
units::UnitsSet the_units):
m_data(the_database),
m_units_set(the_units),
m_dt(HUGE_VAL)
{
}
//! \brief Initialize the stagger at the beginning of the computation
void initialize(VariablesBasePtr var)
{
SaturatedVariablesPtr true_var = cast_var_from_base(var);
true_var->upscaling_variables().setZero();
for (index_t node=0; node<true_var->nb_nodes(); ++node)
{
upscaling_one_node(node, true_var);
}
}
//! \brief Initialize the stagger at the beginning of an iteration
void initialize_timestep(scalar_t dt, VariablesBasePtr var)
{
SaturatedVariablesPtr true_var = cast_var_from_base(var);
m_dt = dt;
for (index_t node=1; node<true_var->nb_nodes(); ++node)
{
upscaling_one_node(node, true_var);
true_var->vel_porosity(node) = 0.0;
}
}
//! \brief Solve the equation for the timestep
StaggerReturnCode restart_timestep(VariablesBasePtr var)
{
SaturatedVariablesPtr true_var = cast_var_from_base(var);
for (index_t node=1; node<true_var->nb_nodes(); ++node)
{
upscaling_one_node(node, true_var);
}
return StaggerReturnCode::ResidualMinimized;
}
scalar_t coeff_diffusion_law(scalar_t porosity)
{
return factor*de_0*((porosity - porosity_r)/(0.4 - porosity_r));
}
void upscaling_one_node(index_t node, SaturatedVariablesPtr true_var)
{
AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node),
m_data, m_units_set);
scalar_t porosity = extractor.porosity();
true_var->vel_porosity(node) += (porosity - true_var->porosity(node))/m_dt;
true_var->porosity(node) = porosity;
true_var->diffusion_coefficient(node) = coeff_diffusion_law(porosity);
}
private:
RawDatabasePtr m_data;
units::UnitsSet m_units_set;
index_t m_id_cshj;
scalar_t m_initial_sat_cshj;
scalar_t m_dt;
};
int main()
{
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
RawDatabasePtr the_database = get_database("../data/cemdata_specmicp.js");
for (auto label: the_database->labels_basis)
std::cout << label << std::endl;
AdimensionalSystemSolverOptions specmicp_options = get_specmicp_options();
units::UnitsSet units_set = specmicp_options.units_set;
AdimensionalSystemConstraints cement_constraints = get_cement_initial_constraints(the_database, units_set);
for (auto label: the_database->labels_minerals)
std::cout << label << std::endl;
AdimensionalSystemSolution cement_solution = get_cement_initial_compo(
the_database,
cement_constraints,
specmicp_options);
AdimensionalSystemSolutionExtractor extractor(cement_solution, the_database, units_set);
std::cout << "Porosity : " << extractor.porosity() << std::endl;
std::cout << "Water volume fraction " << extractor.volume_fraction_water() << std::endl;
AdimensionalSystemConstraints bc_constraints = get_bc_initial_constraints(
the_database,
specmicp_options.units_set
);
AdimensionalSystemSolution bc_solution = get_bc_initial_compo(
the_database,
bc_constraints,
specmicp_options
);
AdimensionalSystemSolutionExtractor extractor2(bc_solution, the_database, units_set);
std::cout << "pH : " << extractor2.pH() << std::endl;
mesh::Mesh1DPtr the_mesh = get_mesh(0.0075);
index_t nb_nodes = the_mesh->nb_nodes();
std::vector<specmicp::AdimensionalSystemConstraints> list_constraints = {cement_constraints,
bc_constraints};
std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition = {cement_solution, bc_solution};
std::vector<int> index_initial_state(nb_nodes, 0);
std::vector<int> index_constraints(nb_nodes, 0);
index_initial_state[0] = 1;
index_constraints[0] = 1 ;
std::vector<index_t> list_fixed_nodes = {0};
satdiff::SaturatedVariablesPtr variables =
satdiff::init_variables(the_mesh, the_database, units_set,
list_fixed_nodes,
list_initial_composition,
index_initial_state);
std::shared_ptr<EquilibriumStagger> equilibrium_stagger =
std::make_shared<satdiff::EquilibriumStagger>(list_constraints, index_constraints, specmicp_options);
std::shared_ptr<SaturatedTransportStagger> transport_stagger =
std::make_shared<satdiff::SaturatedTransportStagger>(variables, list_fixed_nodes);
set_transport_options(transport_stagger->options_solver());
UpscalingStaggerPtr upscaling_stagger =
std::make_shared<CarboUspcaling>(the_database, units_set);
ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger);
transport_stagger->initialize(variables);
equilibrium_stagger->initialize(variables);
upscaling_stagger->initialize(variables);
std::cout << variables->upscaling_variables() << std::endl;
set_reactive_solver_options(solver.get_options());
scalar_t dt = 1.0;
index_t cnt = 0;
Timestepper timestepper(1e-2, 10.0, 30*24*3600, 2);
database::Database dbhandler(the_database);
index_t id_oh = dbhandler.component_label_to_id("HO[-]");
index_t id_hco3 = dbhandler.component_label_to_id("HCO3[-]");
index_t id_ca = dbhandler.component_label_to_id("Ca[2+]");
std::string prefix = "out_carbo_";
std::string suffix = ".dat";
std::ofstream out_c_oh(prefix+"ph"+suffix);
std::ofstream out_c_hco3(prefix+"c_hco3"+suffix);
std::ofstream out_s_ca(prefix+"s_ca"+suffix);
std::ofstream out_s_co2(prefix+"s_co2"+suffix);
std::ofstream out_iter(prefix+"iter"+suffix);
io::print_reactmicp_performance_long_header(&out_iter);
// std::vector<std::ofstream> out_minerals(the_database->nb_mineral);
// for (index_t mineral: the_database->range_mineral())
// {
// const std::string path = prefix+"min_"+the_database->labels_minerals[mineral]+suffix;
// out_minerals[mineral].open(path);
// }
while (timestepper.get_total() < timestepper.get_total_target())
{
Timer step_timer;
step_timer.start();
reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables);
step_timer.stop();
io::print_reactmicp_performance_long(&out_iter, cnt, timestepper.get_total()+dt, solver.get_perfs());
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
{
dt = 0.01;
variables->reset_main_variables();
retcode = solver.solve_timestep(dt, variables);
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
io::print_reactmicp_performance_long(&out_iter, cnt, timestepper.get_total(), solver.get_perfs());
}
if (cnt % 50 == 0)
{
out_c_oh << timestepper.get_total();
out_c_hco3 << timestepper.get_total();
out_s_ca << timestepper.get_total();
out_s_co2 << timestepper.get_total();
//for (index_t mineral: the_database->range_mineral())
// out_minerals[mineral] << timestepper.get_total();
for (index_t node: the_mesh->range_nodes())
{
AdimensionalSystemSolutionExtractor extractor(variables->equilibrium_solution(node), the_database, units_set);
out_c_oh << "\t" << extractor.pH();
out_c_hco3 << "\t" << variables->aqueous_concentration(node, id_hco3, variables->displacement());
out_s_ca << "\t" << variables->solid_concentration(node, id_ca, variables->displacement());
out_s_co2 << "\t" << variables->solid_concentration(node, id_hco3, variables->displacement());
//for (index_t mineral: the_database->range_mineral())
// out_minerals[mineral] << "\t" << extractor.volume_fraction_mineral(mineral);
}
out_c_hco3 << std::endl;
out_c_oh << std::endl;
out_s_ca << std::endl;
out_s_co2 << std::endl;
//for (index_t mineral: the_database->range_mineral())
// out_minerals[mineral] << std::endl;
}
++cnt;
}
print_minerals_profile(the_database, variables, the_mesh, prefix+"profiles_end"+suffix);
io::print_reactmicp_timer(&out_iter, solver.get_timer());
}
AdimensionalSystemSolverOptions get_specmicp_options()
{
AdimensionalSystemSolverOptions options;
options.solver_options.steptol = 1e-10;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.maxstep = 10.0;
options.solver_options.fvectol = 1e-10;
options.solver_options.steptol = 1e-16;
options.solver_options.condition_limit = -1;
options.solver_options.factor_descent_condition = -1;
options.solver_options.threshold_cycling_linesearch = 1e-6;
options.solver_options.non_monotone_linesearch = true;
options.solver_options.use_scaling = true;
options.system_options.non_ideality_tolerance = 1e-14;
options.units_set.length = units::LengthUnit::centimeter;
return options;
}
void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options)
{
transport_options.maximum_iterations = 450;
transport_options.quasi_newton = 3;
transport_options.residuals_tolerance = 1e-6;
transport_options.step_tolerance = 1e-16;
transport_options.absolute_tolerance = 1e-18;
transport_options.alpha = 1.0;
//transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
transport_options.sparse_solver = SparseSolver::SparseLU;
//transport_options.sparse_solver = SparseSolver::GMRES;
transport_options.threshold_stationary_point = 1e-4;
}
void set_reactive_solver_options(reactmicp::solver::ReactiveTransportOptions& reactive_options)
{
reactive_options.maximum_iterations = 50; //500;
reactive_options.residuals_tolerance = 1e-4;
reactive_options.good_enough_tolerance = 1.0;
reactive_options.absolute_residuals_tolerance = 1e-6;
reactive_options.implicit_upscaling = false;
}
RawDatabasePtr get_database(const std::string& path_to_database) {
database::Database the_database(path_to_database);
std::map<std::string, std::string> swapping({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"Al[3+]","Al(OH)4[-]"}
});
the_database.swap_components(swapping);
return the_database.get_database();
}
void print_minerals_profile(
RawDatabasePtr the_database,
SaturatedVariablesPtr variables,
mesh::Mesh1DPtr the_mesh,
const std::string& filepath
)
{
std::ofstream ofile(filepath);
ofile << "Radius";
for (index_t mineral: the_database->range_mineral())
{
ofile << "\t" << the_database->labels_minerals[mineral];
}
ofile << "\t" << "Porosity";
ofile << "\t" << "pH";
ofile << std::endl;
for (index_t node: the_mesh->range_nodes())
{
ofile << the_mesh->get_position(node);
AdimensionalSystemSolutionExtractor extractor(variables->equilibrium_solution(node), the_database, units::UnitsSet());
for (index_t mineral: the_database->range_mineral())
{
ofile << "\t" << extractor.volume_fraction_mineral(mineral);
}
ofile << "\t" << variables->porosity(node);
ofile << "\t" << extractor.pH();
ofile << std::endl;
}
ofile.close();
}
AdimensionalSystemConstraints get_cement_initial_constraints(
RawDatabasePtr the_database,
const units::UnitsSet& units_set)
{
Vector init_min_publi(5);
init_min_publi << 12.9, 1.62, 0.605, 0.131, 0.001;
Vector init_min;
init_min = initial_compo(init_min_publi);
std::cout << "Init min \n ---- \n" << init_min << "\n ---- \n" << std::endl;
scalar_t mult = 1e-6;
//scalar_t mult = 6.73153e-4;
//scalar_t mult = 0.001296315;
init_min *= 6.71146e-4;
Formulation formulation;
formulation.mass_solution = 1.0;
formulation.amount_minerals = {
{"Portlandite", init_min(0)},
{"CSH,jennite", init_min(1)},
{"Monosulfoaluminate", init_min(2)},
{"Ettringite", init_min(3)},
{"Calcite", init_min(4)}
};
formulation.concentration_aqueous =
{
{"Na[+]", mult*1.0298},
{"K[+]", mult*0.0801},
{"Cl[-]", mult*0.0001},
{"HO[-]", mult*1.8298}
},
formulation.minerals_to_keep = {
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2,am",
"Calcite",
"Al(OH)3,am",
"Monosulfoaluminate",
"Tricarboaluminate",
"Monocarboaluminate",
"Hemicarboaluminate",
"Straetlingite",
"Gypsum",
"Ettringite",
"Thaumasite"
};
Dissolver dissolver(the_database);
dissolver.set_units(units_set);
AdimensionalSystemConstraints constraints;
//constraints.set_charge_keeper(1);
constraints.total_concentrations = dissolver.dissolve(formulation);
std::cout << constraints.total_concentrations << std::endl;
constraints.set_saturated_system();
return constraints;
}
AdimensionalSystemConstraints get_bc_initial_constraints(
RawDatabasePtr the_database,
const units::UnitsSet& units_set)
{
const scalar_t rho_w = laws::density_water(units::celsius(25.0), units_set.length, units_set.mass);
const scalar_t nacl = 0.3*rho_w;
const scalar_t kcl = 0.28*rho_w;
database::Database dbhandler(the_database);
Vector total_concentrations(the_database->nb_component);
total_concentrations.setZero();
total_concentrations(dbhandler.component_label_to_id("K[+]")) = kcl;
total_concentrations(dbhandler.component_label_to_id("Cl[-]")) = nacl + kcl;
total_concentrations(dbhandler.component_label_to_id("Na[+]")) = nacl;
AdimensionalSystemConstraints constraints;
constraints.add_fixed_activity_component(dbhandler.component_label_to_id("HO[-]"), -10.3);
constraints.set_charge_keeper(dbhandler.component_label_to_id("HCO3[-]"));
constraints.total_concentrations = total_concentrations;
constraints.set_saturated_system();
return constraints;
}
AdimensionalSystemSolution get_bc_initial_compo(
RawDatabasePtr the_database,
const AdimensionalSystemConstraints& constraints,
const AdimensionalSystemSolverOptions& options
)
{
Vector variables;
AdimensionalSystemSolver solver(the_database,
constraints,
options);
micpsolver::MiCPPerformance perf = solver.solve(variables, true);
if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet)
throw std::runtime_error("Failed to solve initial composition");
std::cout << variables << std::endl;
return solver.get_raw_solution(variables);
}
AdimensionalSystemSolution get_cement_initial_compo(
RawDatabasePtr the_database,
const AdimensionalSystemConstraints& constraints,
const AdimensionalSystemSolverOptions& options
)
{
AdimensionalSystemSolver solver(the_database,
constraints,
options);
Vector variables;
solver.initialise_variables(
variables,
0.4,
{
{"HO[-]", -1.2},
{"Al(OH)4[-]", -4},
{"HCO3[-]", -5},
{"Ca[2+]", -1.6}
});
micpsolver::MiCPPerformance perf = solver.solve(variables);
if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet)
throw std::runtime_error("Failed to solve initial composition");
std::cout << variables << std::endl;
return solver.get_raw_solution(variables);
}
Vector initial_compo(const Vector& publi)
{
Matrix publi_stochio(5,5);
publi_stochio << 1, 0, 0, 0, 0,
9, 6, 0, 0, 0,
3, 0, 1, 1, 0,
3, 0, 1, 3, 0,
0, 0, 0, 0, 1;
Matrix cemdata_stochio(5, 5);
cemdata_stochio << 1, 0, 0, 0, 0,
1.0/0.6, 1.0, 0, 0, 0,
3, 0, 1, 1, 0,
3, 0, 1, 3, 0,
0, 0, 0, 0, 1;
Vector oxyde;
Eigen::ColPivHouseholderQR<Matrix> solver(publi_stochio);
oxyde = solver.solve(publi);
Vector for_cemdata = cemdata_stochio*oxyde;
return for_cemdata;
}
mesh::Mesh1DPtr get_mesh(scalar_t dx)
{
const scalar_t total_length = 7.5/2.0; //cm
const scalar_t height = 20;
index_t nb_nodes = total_length/dx + 1;
return mesh::uniform_axisymmetric_mesh1d(nb_nodes, total_length, dx, height);
}
Event Timeline
Log In to Comment