diff --git a/examples/reactmicp/saturated_react/carbonation.cpp b/examples/reactmicp/saturated_react/carbonation.cpp
index 57de2f4..af2a98c 100644
--- a/examples/reactmicp/saturated_react/carbonation.cpp
+++ b/examples/reactmicp/saturated_react/carbonation.cpp
@@ -1,599 +1,678 @@
#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);
+mesh::Mesh1DPtr get_mesh();
+
+void output_mesh(mesh::Mesh1DPtr the_mesh);
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));
+ if (porosity > 0.7)
+ return factor*4.0*1e-10;
+ else
+ {
+ const scalar_t ratio = ((porosity - porosity_r)/(0.4 - porosity_r));
+ return factor*de_0*std::pow(ratio, 3.32);
+ }
}
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);
+ mesh::Mesh1DPtr the_mesh = get_mesh(0.0025);
+ //mesh::Mesh1DPtr the_mesh = get_mesh(); //(0.0005);
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;
+ output_mesh(the_mesh);
+
+ scalar_t dt = 0.1;
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_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+]");
+ index_t id_calcite = dbhandler.mineral_label_to_id("Calcite");
+
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);
+ std::ofstream out_calcite(prefix+"calcite"+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)
+ if (cnt % 360 == 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();
+ out_calcite << 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());
+ out_calcite << "\t" << extractor.volume_fraction_mineral(id_calcite);
//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;
+ out_calcite << 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"
+ "Ettringite"
};
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 total_length = 0.75/2.0 + dx; //cm
const scalar_t height = 20;
- index_t nb_nodes = total_length/dx + 1;
+ index_t nb_nodes = 151; //total_length/dx + 1;
return mesh::uniform_axisymmetric_mesh1d(nb_nodes, total_length, dx, height);
}
+
+mesh::Mesh1DPtr get_mesh()
+{
+ Vector radius(40);
+
+ const scalar_t height = 20;
+
+ radius <<
+ 3.751,
+ 3.750,
+ 3.749,
+ 3.748,
+ 3.747,
+ 3.745,
+ 3.740,
+ 3.735,
+ 3.730,
+ 3.720,
+ 3.710,
+ 3.700,
+ 3.675,
+ 3.650,
+ 3.625,
+ 3.600,
+ 3.575,
+ 3.550,
+ 3.525,
+ 3.500,
+ 3.475,
+ 3.450,
+ 3.425,
+ 3.400,
+ 3.375,
+ 3.350,
+ 3.325,
+ 3.300,
+ 3.250,
+ 3.200,
+ 3.150,
+ 3.100,
+ 3.050,
+ 3.000,
+ 2.950,
+ 2.900,
+ 2.850,
+ 2.800,
+ 2.750,
+ 2.700;
+
+
+ return mesh::axisymmetric_mesh1d(radius, height);
+}
+
+void output_mesh(mesh::Mesh1DPtr the_mesh)
+{
+ std::ofstream output_mesh;
+ output_mesh.open("out_carbo_mesh.dat");
+
+ output_mesh << "Node\tPosition" << std::endl;
+ for (index_t node: the_mesh->range_nodes())
+ {
+ output_mesh << node << "\t" << the_mesh->get_position(node) << std::endl;
+ }
+ output_mesh.close();
+
+}