Homec4science
Diffusion SpecMiCP / ReactMiCP (configuration) Browse

specmicp/examples/reactmicp/saturated_react/react_leaching.cpp
532031755403configuration

Last Change
Raw File

react_leaching.cpp
View Options

#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp"
#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
#include "reactmicp/solver/reactive_transport_solver.hpp"
#include "reactmicp/systems/saturated_react/variables.hpp"
#include "reactmicp/systems/saturated_react/transport_stagger.hpp"
#include "reactmicp/systems/saturated_react/equilibrium_stagger.hpp"
#include "reactmicp/systems/saturated_react/init_variables.hpp"
#include "dfpm/meshes/axisymmetric_uniform_mesh1d.hpp"
#include "dfpm/meshes/axisymmetric_mesh1d.hpp"
#include "utils/log.hpp"
#include <iostream>
#include <fstream>
#include "utils/io/reactive_transport.hpp"
using namespace specmicp;
// Initial conditions
// ==================
specmicp::RawDatabasePtr leaching_database()
{
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.remove_gas_phases();
thedatabase.swap_components(swapping);
return thedatabase.get_database();
}
AdimensionalSystemSolution initial_leaching_pb(
RawDatabasePtr raw_data,
scalar_t mult,
units::UnitsSet the_units)
{
Formulation formulation;
scalar_t m_c3s = mult*0.6;
scalar_t m_c2s = mult*0.2;
scalar_t m_c3a = mult*0.1;
scalar_t m_gypsum = mult*0.1;
scalar_t wc = 0.5;
scalar_t m_water = wc*1.0e-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}
};
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
constraints.charge_keeper = 1;
constraints.set_saturated_system();
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 100.0;
options.solver_options.max_iter = 200;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = true;
options.solver_options.factor_descent_condition = -1.0;
options.solver_options.factor_gradient_search_direction = 100;
options.solver_options.projection_min_variable = 1e-9;
options.solver_options.fvectol = 1e-10;
options.solver_options.steptol = 1e-14;
options.solver_options.non_monotone_linesearch = true;
options.system_options.non_ideality_tolerance = 1e-10;
options.units_set = the_units;
specmicp::Vector x;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
solver.initialise_variables(x, 0.8, -3.0);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
std::cout << "Initial return code : " << (int) perf.return_code << std::endl;
if ((int) perf.return_code <= (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet)
{
throw std::runtime_error("Initial solution has not converged");
}
return solver.get_raw_solution(x);
}
AdimensionalSystemSolution initial_blank_leaching_pb(
RawDatabasePtr raw_data,
scalar_t mult,
units::UnitsSet the_units)
{
Formulation formulation;
scalar_t m_c3s = mult*0.6;
scalar_t m_c2s = mult*0.2;
scalar_t m_c3a = mult*0.1;
scalar_t m_gypsum = mult*0.1;
scalar_t wc = 0.8;
scalar_t m_water = wc*1.0e-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 = {
{"SiO2,am", mult},
};
// formulation.concentration_aqueous = {
// {"Ca[2+]", 1e-7},
// {"Al(OH)4[-]", 1e-7},
// {"SO4[2-]", 1e-7}
// };
formulation.extra_components_to_keep = {"Ca[2+]", "SO4[2-]", "Al(OH)4[-]"};
specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
constraints.charge_keeper = 1;
constraints.set_saturated_system();
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 10.0;
options.solver_options.max_iter = 200;
options.solver_options.maxiter_maxstep = 200;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = true;
options.solver_options.factor_descent_condition = -1.0;
options.solver_options.factor_gradient_search_direction = 100;
options.solver_options.projection_min_variable = 1e-9;
options.solver_options.fvectol = 1e-4;
options.solver_options.steptol = 1e-14;
options.solver_options.non_monotone_linesearch = true;
options.system_options.non_ideality_tolerance = 1e-10;
options.units_set = the_units;
Vector x;
AdimensionalSystemSolver solver(raw_data, constraints, options);
solver.initialise_variables(x, 0.75, -4.0);
micpsolver::MiCPPerformance perf = solver.solve(x, false);
std::cout << "Initial return code : " << (int) perf.return_code << std::endl;
if ((int) perf.return_code <= (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet)
{
throw std::runtime_error("Initial solution has not converged");
}
AdimensionalSystemSolution solution = solver.get_raw_solution(x);
AdimensionalSystemSolutionModificator extractor(solution, raw_data, the_units);
extractor.remove_solids();
Vector totconc(raw_data->nb_component);
totconc(0) = extractor.density_water()*extractor.total_aqueous_concentration(0);
for (index_t component: raw_data->range_aqueous_component())
{
totconc(component) = 0.01*extractor.density_water()*extractor.total_aqueous_concentration(component);
}
constraints.total_concentrations = totconc;
solver = AdimensionalSystemSolver(raw_data, constraints, options);
perf = solver.solve(x, false);
specmicp_assert((int) perf.return_code > (int) micpsolver::MiCPSolverReturnCode::NotConvergedYet);
return solver.get_raw_solution(x);
}
// Upscaling
// =========
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::solver;
using namespace specmicp::reactmicp::systems::satdiff;
class LeachingUspcaling: public solver::UpscalingStaggerBase
{
public:
LeachingUspcaling(RawDatabasePtr the_database,
units::UnitsSet the_units):
m_data(the_database),
m_units_set(the_units)
{
}
//! \brief Initialize the stagger at the beginning of the computation
void initialize(VariablesBasePtr var)
{
SaturatedVariablesPtr true_var = cast_var_from_base(var);
database::Database db_handler(m_data);
m_dt = 1;
m_id_cshj = db_handler.mineral_label_to_id("CSH,j");
m_initial_sat_cshj = true_var->equilibrium_solution(1).main_variables(m_data->nb_component+m_id_cshj);
//m_initial_sat_cshj = 0.36;
std::cout << m_initial_sat_cshj << std::endl;
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;
}
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->equilibrium_solution(node).main_variables(m_data->nb_component)/m_initial_sat_cshj*0.05;
//std::cout << "porosity : " << porosity
// << " - saturation water" << true_var->equilibrium_solution(node).main_variables(0) << std::endl;
true_var->vel_porosity(node) += (porosity - true_var->porosity(node))/m_dt;
true_var->porosity(node) = porosity;
true_var->diffusion_coefficient(node) = std::min(1e4*std::exp(9.85*porosity-29.08),1e-6);
}
private:
RawDatabasePtr m_data;
units::UnitsSet m_units_set;
index_t m_id_cshj;
scalar_t m_initial_sat_cshj;
scalar_t m_dt;
};
//
// Main
// ====
int main()
{
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
RawDatabasePtr raw_data = leaching_database();
units::UnitsSet the_units;
the_units.length = units::LengthUnit::centimeter;
AdimensionalSystemSolution initial = initial_leaching_pb(raw_data, 0.006, the_units);
AdimensionalSystemSolutionModificator extractor(initial, raw_data, the_units);
Database dbhandler(raw_data);
index_t id_ca = dbhandler.component_label_to_id("Ca[2+]");
index_t id_si = dbhandler.component_label_to_id("SiO(OH)3[-]");
std::cout << extractor.total_solid_concentration(id_ca) << " - " << extractor.porosity() << std::endl;
extractor.scale_total_concentration(id_ca, 0.015);
std::cout << extractor.total_solid_concentration(id_ca) << " - " << extractor.porosity() << std::endl;
std::cout << initial.main_variables << std::endl;
AdimensionalSystemSolution blank = initial_blank_leaching_pb(raw_data, 0.005, the_units);
std::cout << blank.main_variables << std::endl;
scalar_t radius = 3.5; //cm
scalar_t height = 8.0; //cm
scalar_t dx = 0.0025;
index_t additional_nodes = 1;
radius = radius + additional_nodes*dx;
index_t nb_nodes = 100;
// mesh::Mesh1DPtr the_mesh = mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, dx, height);
mesh::Mesh1DPtr the_mesh = mesh::uniform_axisymmetric_mesh1d(nb_nodes, radius, dx, height);
std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition = {initial, blank};
std::vector<int> index_initial_state(nb_nodes, 0);
for (int i=0; i< additional_nodes; ++i)
index_initial_state[i] = 1;
std::vector<specmicp::index_t> list_fixed_nodes = {0, };
systems::satdiff::SaturatedVariablesPtr variables =
systems::satdiff::init_variables(the_mesh, raw_data, the_units,
list_fixed_nodes,
list_initial_composition,
index_initial_state);
specmicp::AdimensionalSystemConstraints constraints;
constraints.charge_keeper = 1;
constraints.set_saturated_system();
constraints.inert_volume_fraction = 0.0;
AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 100.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-16;
options.solver_options.fvectol = 1e-8;
options.solver_options.steptol = 1e-12;
options.solver_options.non_monotone_linesearch = true;
options.system_options.non_ideality_tolerance = 1e-12;
options.system_options.non_ideality_max_iter = 100;
options.units_set = the_units;
ChemistryStaggerPtr equilibrium_stagger =
std::make_shared<reactmicp::systems::satdiff::EquilibriumStagger>(nb_nodes, constraints, options);
std::shared_ptr<reactmicp::systems::satdiff::SaturatedTransportStagger> transport_stagger =
std::make_shared<reactmicp::systems::satdiff::SaturatedTransportStagger>(variables, list_fixed_nodes);
dfpmsolver::ParabolicDriverOptions& transport_options = transport_stagger->options_solver();
transport_options.maximum_iterations = 450;
transport_options.quasi_newton = 3;
transport_options.residuals_tolerance = 1e-8;
transport_options.step_tolerance = 1e-10;
transport_options.alpha = 1;
transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
transport_options.sparse_solver = SparseSolver::SparseLU;
//transport_options.sparse_solver = SparseSolver::GMRES;
transport_options.threshold_stationary_point = 1e-2;
UpscalingStaggerPtr upscaling_stagger =
std::make_shared<LeachingUspcaling>(raw_data, the_units);
ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger);
transport_stagger->initialize(variables);
equilibrium_stagger->initialize(variables);
upscaling_stagger->initialize(variables);
solver.get_options().maximum_iterations = 500;
solver.get_options().residuals_tolerance = 1e-2;
solver.get_options().absolute_residuals_tolerance = 1e-18;
solver.get_options().implicit_upscaling = true;
std::ofstream output_iter("out_leaching_iter.dat");
std::ofstream output_c_ca, output_s_ca, output_s_si;
std::ofstream output_poro;
std::ofstream output_loss_ca;
output_c_ca.open("out_c_ca.dat");
output_s_ca.open("out_s_ca.dat");
output_s_si.open("out_s_si.dat");
output_poro.open("out_poro.dat");
output_loss_ca.open("out_loss_ca.dat");
scalar_t initial_ca = 0;
for (index_t node: the_mesh->range_nodes())
{
initial_ca += variables->solid_concentration(node, id_ca, variables->displacement())
* the_mesh->get_volume_cell(node);
}
scalar_t total = 0.0;
scalar_t dt=2.0;
io::print_reactmicp_performance_long_header(&output_iter);
int cnt =0;
while (total < 25*3600*24)
{
solver.solve_timestep(dt, variables);
total += dt;
io::print_reactmicp_performance_long(&output_iter, cnt, total, solver.get_perfs());
std::cout << "time : " << total/3600/24 << std::endl
<< "Iter : " << solver.get_perfs().nb_iterations << std::endl
<< "Residuals : " << solver.get_perfs().residuals << std::endl;
if (cnt % 5000 == 0)
{
scalar_t time = total/3600/24;
scalar_t sqrt_time = std::sqrt(time);
output_c_ca << time;
output_s_ca << time;
output_s_si << time;
output_loss_ca << time << "\t" << sqrt_time;
output_poro << time;
scalar_t amount_ca = 0.0;
for (index_t node: the_mesh->range_nodes())
{
amount_ca += variables->solid_concentration(node, id_ca, variables->displacement())
* the_mesh->get_volume_cell(node);
output_c_ca << "\t" << variables->aqueous_concentration(node, id_ca, variables->displacement());
output_s_ca << "\t" << variables->solid_concentration( node, id_ca, variables->displacement());
output_s_si << "\t" << variables->solid_concentration( node, id_si, variables->displacement());
output_poro << "\t" << variables->porosity(node);
}
output_poro << std::endl;
output_loss_ca << "\t" << (initial_ca - amount_ca)/1.75929e-2 << std::endl;
output_c_ca << std::endl;
output_s_ca << std::endl;
output_s_si << std::endl;
}
++cnt;
}
io::print_reactmicp_timer(&std::cout, solver.get_timer());
}
c4science · Help