Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F93730444
react_leaching.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
Sun, Dec 1, 01:11
Size
16 KB
Mime Type
text/x-c++
Expires
Tue, Dec 3, 01:11 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
22695236
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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>
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(raw_data->nb_component+raw_data->nb_mineral);
x.setZero();
x(0) = 0.8;
x.segment(1, raw_data->nb_component-1).setConstant(-3.0);
specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
specmicp::micpsolver::MiCPPerformance perf = solver.solve(x, false);
specmicp_assert((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
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(raw_data->nb_component+raw_data->nb_mineral);
x.setZero();
x(0) = 0.75;
x.segment(1, raw_data->nb_component-1).setConstant(-4.0);
AdimensionalSystemSolver solver(raw_data, constraints, options);
micpsolver::MiCPPerformance perf = solver.solve(x, false);
specmicp_assert((int) perf.return_code > (int) micpsolver::MiCPSolverReturnCode::NotConvergedYet);
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;
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 = 50.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-10;
options.solver_options.fvectol = 1e-10;
options.solver_options.steptol = 1e-14;
options.solver_options.non_monotone_linesearch = false;
options.system_options.non_ideality_tolerance = 1e-14;
options.system_options.non_ideality_max_iter = 100;
options.units_set = the_units;
ChemistryStaggerPtr equilibrium_stagger =
std::make_shared<reactmicp::systems::satdiff::EquilibriumStagger>(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_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=1.0;
int cnt =0;
while (total < 25*3600*24)
{
solver.solve_timestep(dt, variables);
//std::cout<< variables->displacement() << std::endl;
total += dt;
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;
}
}
Event Timeline
Log In to Comment