Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/examples/reactmicp/equilibrium_curve.cpp
edda4626caa0reactivetransport

Last Change
Raw File

equilibrium_curve.cpp
View Options

#include <iostream>
#include "utils/log.hpp"
#include "reactmicp/equilibrium_curve/chemistry.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "reactmicp/equilibrium_curve/eqcurve_extractor.hpp"
#include "reactmicp/equilibrium_curve/eqcurve_coupler.hpp"
#include "dfpm/meshes/axisymmetric_uniform_mesh1d.hpp"
#include "dfpm/meshes/uniform_mesh1d.hpp"
#include "reactmicp/equilibrium_curve/eqcurve_solid_transport.hpp"
specmicp::Matrix test_chemistry() {
specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
});
thedatabase.swap_components(swapping);
thedatabase.remove_gas_phases();
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::Formulation formulation;
specmicp::scalar_t mult = 6.5e3;
specmicp::scalar_t m_c3s = mult*0.7;
specmicp::scalar_t m_c2s = mult*0.3;
specmicp::scalar_t wc = 0.5;
specmicp::scalar_t m_water = wc*1e-3*(
m_c3s*(3*56.08+60.08)
+ m_c2s*(2*56.06+60.08)
);
formulation.mass_solution = m_water;
formulation.amount_minerals = {
{"C3S", m_c3s},
{"C2S", m_c2s},
};
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
specmicp::index_t id_h2o = thedatabase.component_label_to_id("H2O");
specmicp::index_t id_ho = thedatabase.component_label_to_id("HO[-]");
specmicp::index_t id_ca = thedatabase.component_label_to_id("Ca[2+]");
specmicp::AdimensionalSystemBC conditions(total_concentrations);
conditions.charge_keeper = id_ho;
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 10.0;
options.solver_options.max_iter = 100;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = false;
options.solver_options.factor_descent_condition = -1;
options.solver_options.factor_gradient_search_direction = 100;
options.solver_options.projection_min_variable = 1e-9;
options.solver_options.fvectol = 1e-6;
options.solver_options.steptol = 1e-14;
options.system_options.non_ideality_tolerance = 1e-10;
specmicp::reactmicp::eqcurve::EquilibriumCurveSpeciation spec_solver(raw_data, conditions, id_ca, options);
return spec_solver.get_equilibrium_curve(0.05, -500);
}
specmicp::Matrix test_chemistry_with_al()
{
specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
std::map<std::string, std::string> swapping ({
{"H[+]","HO[-]"},
{"Si(OH)4", "SiO(OH)3[-]"},
{"Al[3+]","Al(OH)4[-]"}
});
thedatabase.swap_components(swapping);
thedatabase.remove_gas_phases();
specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
specmicp::Formulation formulation;
specmicp::scalar_t mult = 6.5e3;
specmicp::scalar_t m_c3s = mult*0.6;
specmicp::scalar_t m_c2s = mult*0.2;
specmicp::scalar_t m_c3a = mult*0.10;
specmicp::scalar_t m_gypsum = mult*0.10;
specmicp::scalar_t wc = 0.8;
specmicp::scalar_t m_water = wc*1e-3*(
m_c3s*(3*56.08+60.08)
+ m_c2s*(2*56.06+60.08)
+ m_c3a*(3*56.08+101.96)
+ m_gypsum*(56.08+80.06+2*18.02)
);
formulation.mass_solution = m_water;
formulation.amount_minerals = {
{"C3S", m_c3s},
{"C2S", m_c2s},
{"C3A", m_c3a},
{"Gypsum", m_gypsum}
};
formulation.minerals_to_keep = {
"Portlandite",
"CSH,jennite",
"CSH,tobermorite",
"SiO2,am",
"Al(OH)3,am",
"Monosulfoaluminate",
"Straetlingite",
"Gypsum",
"Ettringite",
};
for (specmicp::index_t component: raw_data->range_component())
{
std::cout << raw_data->labels_basis[component] << std::endl;
}
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
//specmicp::index_t id_h2o = thedatabase.component_label_to_id("H2O");
specmicp::index_t id_ho = thedatabase.component_label_to_id("HO[-]");
specmicp::index_t id_ca = thedatabase.component_label_to_id("Ca[2+]");
specmicp::AdimensionalSystemBC conditions(total_concentrations);
conditions.charge_keeper = id_ho;
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 20.0;
options.solver_options.max_iter = 100;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = false;
options.solver_options.factor_descent_condition = -1;
options.solver_options.factor_gradient_search_direction = 100;
options.solver_options.projection_min_variable = 1e-9;
options.solver_options.fvectol = 1e-6;
options.solver_options.steptol = 1e-14;
options.system_options.non_ideality_tolerance = 1e-10;
specmicp::reactmicp::eqcurve::EquilibriumCurveSpeciation spec_solver(raw_data, conditions, id_ca, options);
return spec_solver.get_equilibrium_curve(0.05, -500.0);
}
void test_eqcurve_extractor()
{
specmicp::reactmicp::eqcurve::EquilibriumCurveExtractor extract(test_chemistry_with_al());
specmicp::index_t index = extract.find_point(111.0);
std::cout << "111.0 \t " << extract.totsolid_concentration(index) << "\t"
<< extract.totaq_concentration(index) << "\t"
<< extract.porosity(index) << "\t"
<< extract.diffusion_coefficient(index) << std::endl;
}
void test_diffeqcurve()
{
specmicp::Matrix eqcurve = test_chemistry_with_al();
eqcurve.col(0) *= 1e-6; //mol/m3 -> mol/cm3
eqcurve.col(1) *= 1e-3; //mol/kg -> mol/cm3
std::cout << eqcurve << std::endl;
specmicp::scalar_t radius = 3.5; //cm
specmicp::scalar_t height = 8.0; //cm
specmicp::scalar_t dx = 0.0025;
specmicp::index_t additional_nodes = 1;
radius = radius + additional_nodes*dx;
specmicp::index_t nb_nodes =10+additional_nodes;
specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, dx, height);
specmicp::dfpmsolver::ParabolicDriverOptions options;
options.step_tolerance = 1e-10;
options.residuals_tolerance = 1e-8;
options.quasi_newton = 1;
specmicp::reactmicp::eqcurve::EquilibriumCurveCoupler solver(eqcurve, the_mesh, options);
specmicp::scalar_t sum_0 =0;
for (specmicp::index_t node=0; node< the_mesh->nb_nodes(); ++node)
{
sum_0 += solver.solid_concentrations()(node)*the_mesh->get_volume_cell(node);
std::cout << the_mesh->get_volume_cell(node) << std::endl;
}
specmicp::scalar_t dt = 0.5;
specmicp::scalar_t total = 0;
std::cout << total << "\t" << 0.0 << "\t" << sum_0 << "\t" << 0.0 << std::endl;
specmicp::index_t k=0;
while (total < 0.4)
{
solver.run_step(dt);
total += dt/3600/24;
if (k % 5000 == 0)
{
specmicp::scalar_t sum =0;
for (specmicp::index_t node=0; node< the_mesh->nb_nodes(); ++node)
{
sum += solver.solid_concentrations()(node)*the_mesh->get_volume_cell(node);
}
std::cout << total << "\t" << std::sqrt(total) << "\t" << sum << "\t" << (sum_0 -sum)/(1.75929e-2) << std::endl;
}
++k;
}
//std::cout << solver.solid_concentrations() << std::endl;
}
void test_eqcurve_solid()
{
//specmicp::Matrix eq_curve = test_chemistry();
specmicp::Matrix eq_curve = test_chemistry_with_al();
eq_curve.col(0) *= 1e-6; //mol/m3 -> mol/cm3
eq_curve.col(1) *= 1e-3; //mol/kg -> mol/cm3
for (specmicp::index_t ind=1; ind<eq_curve.rows(); ++ind)
{
if (eq_curve(ind, 1) >= eq_curve(ind-1, 1))
eq_curve(ind, 1) = eq_curve(ind-1, 1);
}
std::cout << eq_curve << std::endl;
specmicp::scalar_t radius = 3.5; //cm
specmicp::scalar_t height = 8.0; //cm
specmicp::scalar_t dx = 0.005;
specmicp::index_t additional_nodes = 1;
radius = radius + additional_nodes*dx;
specmicp::index_t nb_nodes = 50+additional_nodes;
specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, dx, height);
//specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::uniform_mesh1d(nb_nodes, dx, 5);
specmicp::dfpmsolver::ParabolicDriverOptions options;
options.step_tolerance = 1e-12;
options.residuals_tolerance = 1e-6;
options.sparse_solver = specmicp::SparseSolver::GMRES;
//options.linesearch = specmicp::dfpmsolver::ParabolicLinesearch::Strang;
options.alpha = 1.0;
options.quasi_newton = 1;
options.maximum_step_length = 10;
specmicp::reactmicp::eqcurve::SolidDiffusion program(the_mesh, eq_curve, {0,});
specmicp::dfpmsolver::ParabolicDriver< specmicp::reactmicp::eqcurve::SolidDiffusion> solver(program);
solver.get_options() = options;
solver.set_scaling(specmicp::Vector::Constant(program.get_neq(), 1e6));
specmicp::Vector variables(nb_nodes);
variables(0) = eq_curve(eq_curve.rows()-10, 0);
for (specmicp::index_t node=1; node<the_mesh->nb_nodes(); ++node)
{
variables(node) = eq_curve(10, 0);
}
specmicp::scalar_t sum_0 =0;
for (specmicp::index_t node=0; node< the_mesh->nb_nodes(); ++node)
{
sum_0 += variables(node)*the_mesh->get_volume_cell(node);
std::cout << the_mesh->get_volume_cell(node) << std::endl;
}
specmicp::scalar_t dt = 10.0;
specmicp::scalar_t total = 0;
std::cout << total << "\t" << 0.0 << "\t" << sum_0 << "\t" << 0.0 << std::endl;
specmicp::index_t k=0;
while (total < 50)
{
//std::cout << " ==== TIMESTEP === " << std::endl;
solver.solve_timestep(dt, variables);
for (int node=0; node<the_mesh->nb_nodes(); ++node)
if (variables(node) < 1e-6) variables(node) = 0;
//std::cout << solver.get_perfs().nb_iterations << std::endl;
total += dt/3600/24;
if (k % 5000 == 0)
{
specmicp::scalar_t sum =0;
for (specmicp::index_t node=0; node< the_mesh->nb_nodes(); ++node)
{
sum += variables(node)*the_mesh->get_volume_cell(node);
}
std::cout << total << "\t" << std::sqrt(total) << "\t" << sum << "\t" << (sum_0 -sum)/(1.75929e-2) << std::endl;
}
++k;
}
std::cout << variables << std::endl;
}
void test_interpolator()
{
specmicp::Matrix mat(5,4);
mat << 1, 1, 1, 1,
2, 1, 2, 0,
3, 1, 3, -1,
4, 1, 4, -2,
5, 1, 5, -3;
specmicp::reactmicp::eqcurve::EquilibriumCurveExtractor interpolator(mat);
std::cout << " " << interpolator.slope(0, 1) << " ?== " << 0 << std::endl;
std::cout << " " << interpolator.slope(0, 2) << " ?== " << 1 << std::endl;
std::cout << interpolator.slope(0, 3) << " ?== " << -1 << std::endl;
std::cout << interpolator.find_point(1.5) << " ? == " << 0 << std::endl;
std::cout << interpolator.find_point(3.5) << " ? == " << 2 << std::endl;
std::cout << interpolator.interpolate(2, 3.5, 2) << " ? == " << 3.5 << std::endl;
std::cout << interpolator.interpolate(2, 3.5, 3) << " ? == " << -1.5 << std::endl;
specmicp::Matrix mat2(5,4);
mat2 << 5, 1, 1, 1,
4, 1, 2, 0,
3, 1, 3, -1,
2, 1, 4, -2,
1, 1, 5, -3;
specmicp::reactmicp::eqcurve::EquilibriumCurveExtractor interpolator2(mat2);
std::cout << " " << interpolator2.slope(0, 1) << " ?== " << 0 << std::endl;
std::cout << " " << interpolator2.slope(0, 2) << " ?== " << -1 << std::endl;
std::cout << interpolator2.slope(0, 3) << " ?== " << +1 << std::endl;
std::cout << interpolator2.find_point(1.5) << " ? == " << 3 << std::endl;
std::cout << interpolator2.find_point(3.5) << " ? == " << 1 << std::endl;
std::cout << interpolator2.interpolate(1, 3.5, 2) << " ? == " << 2.5 << std::endl;
std::cout << interpolator2.interpolate(1, 3.5, 3) << " ? == " << -0.5 << std::endl;
std::cout << interpolator2.interpolate(4, 1.0, 2) << " ? == " << 5 << std::endl;
std::cout << interpolator2.interpolate(4, 1.0, 3) << " ? == " << -3 << std::endl;
}
int main()
{
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
specmicp::logger::ErrFile::stream() = &std::cerr;
//test_chemistry();
//std::cout << test_chemistry_with_al() << std::endl;
//test_eqcurve_extractor();
//test_diffeqcurve();
test_eqcurve_solid();
test_interpolator();
return EXIT_SUCCESS;
}
c4science · Help