diff --git a/examples/reactmicp/leaching.cpp b/examples/reactmicp/leaching.cpp index 1a4c8be..83236a8 100644 --- a/examples/reactmicp/leaching.cpp +++ b/examples/reactmicp/leaching.cpp @@ -1,518 +1,561 @@ #include #include "database/database.hpp" #include "specmicp/reaction_path.hpp" #include "reactmicp/mesh.hpp" #include "reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp" #include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp" #include "reactmicp/systems/saturated_diffusion/boundary_conditions.hpp" #include "reactmicp/systems/saturated_diffusion/reactive_transport_neutrality_solver.hpp" #include "reactmicp/systems/saturated_diffusion/transport_neutrality_parameters.hpp" #include #include using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems::siasaturated; const double M_CaO = 56.08; const double M_SiO2 = 60.09; const double M_Al2O3 = 101.96; const double M_SO3 = 80.06; void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species) { Eigen::MatrixXd compo_in_oxyde(4, 4); Vector molar_mass(4); molar_mass << 1.0/M_CaO, 1.0/M_SiO2, 1.0/M_Al2O3, 1.0/M_SO3; compo_oxyde = compo_oxyde.cwiseProduct(molar_mass); //C3S C2S C3A Gypsum compo_in_oxyde << 3, 2, 3, 1, // CaO 1, 1, 0, 0, // SiO2 0, 0, 1, 0, // Al2O3 0, 0, 0, 1; // SO3 Eigen::FullPivLU solver(compo_in_oxyde); compo_species = solver.solve(compo_oxyde); } std::shared_ptr set_database() { specmicp::database::Database database("../data/cemdata_specmicp.js"); std::map swapping ({ {"H[+]","HO[-]"}, {"Si(OH)4", "SiO(OH)3[-]"}, {"Al[3+]","Al(OH)4[-]"} }); database.swap_components(swapping); //std::vector to_keep = {"H2O", "HO[-]", "Ca[2+]", "SiO(OH)3[-]"}; //database.keep_only_components(to_keep); return database.get_database(); } EquilibriumState get_initial_state( std::shared_ptr database, Vector& compo, double mult) { std::shared_ptr model = std::make_shared(); compo = mult*compo; const double m_c3s = compo(0); const double m_c2s = compo(1); const double m_c3a = compo(2); const double m_gypsum = compo(3); const double wc = 1.0; const double m_water = wc*( (3*M_CaO+M_SiO2)*m_c3s + (2*M_CaO+M_SiO2)*m_c2s + (3*M_CaO+M_Al2O3)*m_c3a + (M_CaO + M_SO3)*m_gypsum )*1e-3; model->amount_aqueous = { {"H2O", specmicp::reaction_amount_t(m_water/database->molar_mass_basis_si(0), 0.0)}, {"HO[-]", specmicp::reaction_amount_t(0.0, 0.0)}, }; model->amount_minerals = { {"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)}, {"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)}, {"C3A", specmicp::reaction_amount_t(m_c3a, 0.0)}, {"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0.0)}, }; model->minerals_to_keep = { "Portlandite", "CSH,jennite", "CSH,tobermorite", "SiO2,am", "Al(OH)3,am", "Monosulfoaluminate", "Straetlingite", "Gypsum", "Ettringite" }; Eigen::VectorXd x; specmicp::ReactionPathDriver driver(model, database, x); //driver.get_options().solver_options.penalization_factor =1.0; driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB; //driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min; driver.get_options().solver_options.use_scaling = false; driver.get_options().solver_options.max_factorization_step = 4; driver.get_options().solver_options.factor_gradient_search_direction = 200; driver.get_options().solver_options.maxstep = 50; driver.get_options().solver_options.maxiter_maxstep = 50; driver.get_options().solver_options.max_iter = 100; driver.get_options().allow_restart = true; specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x); if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized) throw std::runtime_error("Error - system cannot be solved"); return driver.get_current_solution(); } EquilibriumState get_blank_state(std::shared_ptr database) { std::shared_ptr model = std::make_shared(); - const double m_c3s = 1e-6; - const double m_c2s = 1e-6; - const double m_c3a = 1e-6; - const double m_gypsum = 1e-6; - const double m_water = 1; + const double m_c3s = 1e-8; + const double m_c2s = 1e-8; + const double m_c3a = 10e-8; + const double m_gypsum = 10e-8; + const double m_water = 1.0; model->amount_aqueous = { {"H2O", specmicp::reaction_amount_t(m_water/database->molar_mass_basis_si(0), 0.0)}, {"HO[-]", specmicp::reaction_amount_t(-0.0, 0.0)}, {"SiO(OH)3[-]", specmicp::reaction_amount_t(0.0, 0.0)} }; model->amount_minerals = { - {"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)}, + //{"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)}, {"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)}, {"C3A", specmicp::reaction_amount_t(m_c3a, 0.0)}, {"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0.0)}, - //{"SiO2,am", specmicp::reaction_amount_t(0.0, 0.0)} + //{"SiO2,am", specmicp::reaction_amount_t(0.0001, 0.0)} }; - Eigen::VectorXd x; + Eigen::VectorXd x(database->nb_component+database->nb_mineral); + //x.setZero(); + //x(0) = 1.0; + //x.segment(1, database->nb_component-1).setConstant(-6); specmicp::ReactionPathDriver driver(model, database, x); //driver.get_options().solver_options.penalization_factor =1.0; driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB; //driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min; driver.get_options().solver_options.use_scaling = false; driver.get_options().solver_options.max_factorization_step = 1; driver.get_options().solver_options.factor_gradient_search_direction = 50; driver.get_options().solver_options.maxstep = 50; driver.get_options().solver_options.maxiter_maxstep = 50; driver.get_options().solver_options.max_iter = 100; driver.get_options().allow_restart = true; specmicp::micpsolver::MiCPPerformance perf = driver.one_step(x); if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized) throw std::runtime_error("Error - blank system cannot be solved"); + + EquilibriumState solution = driver.get_current_solution(); + std::cout << solution.main_variables() << std::endl; + solution.remove_minerals(); + + + + return driver.get_current_solution(); } double get_dt( std::shared_ptr parameters, std::shared_ptr themesh ) { const scalar_t max_d = parameters->diffusion_coefficient(1); scalar_t dt = std::pow(themesh->get_dx(0), 2)/max_d/3; - if (dt > 10) return 10; + if (dt > 5.0) return 5.0; //return 10.0; //return dt; //if (dt<3600*24) dt=3600*24; //if (dt>2*3600*24) dt=2*3600*24; //return 900; } int leaching() { Vector compo_oxyde(4), compo_species(4); compo_oxyde << 62.9, 20.6, 5.8, 3.1; compo_from_oxyde(compo_oxyde, compo_species); std::shared_ptr thedatabase = set_database(); scalar_t radius = 3.5; //cm - scalar_t height = 8; //cm - index_t nb_nodes = 80; + scalar_t height = 8.0; //cm + index_t nb_nodes = 141; mesh::Mesh1DPtr themesh = mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, height); std::shared_ptr parameters = - std::make_shared(nb_nodes, 1e4*std::exp(9.95*0.15-29.08), 0.15); - parameters->diffusion_coefficient(0) = 2.2e-5; + std::make_shared(nb_nodes, 1e4*std::exp(9.95*0.25-29.08), 0.25); + parameters->diffusion_coefficient(0) = 2.2e-5; // cm^2/s parameters->porosity(0) = 1.0; EquilibriumState initial_state = get_initial_state(thedatabase, compo_species, 0.5); SIABoundaryConditions bcs(nb_nodes); bcs.list_initial_states.push_back(initial_state); bcs.list_initial_states.push_back(get_blank_state(thedatabase)); bcs.bs_types[0] = BoundaryConditionType::FixedComposition; bcs.initial_states[0] = 1; database::Database dbhandler(thedatabase); index_t id_ca = dbhandler.component_label_to_id("Ca[2+]"); index_t id_si = dbhandler.component_label_to_id("SiO(OH)3[-]"); index_t id_oh = dbhandler.component_label_to_id("HO[-]"); + for (index_t mineral: thedatabase->range_mineral()) + { + std::cout << thedatabase->labels_minerals[mineral] << std::endl; + } + + std::cout << "Volume first cell " << themesh->get_volume_cell(1) << std::endl; + + //return EXIT_SUCCESS; + std::ofstream output_ca, output_aqca, output_si, output_poro, output_ph; + std::ofstream output_saca, output_sasi; output_ca.open("out_ca.dat"); output_aqca.open("out_aqca.dat"); output_si.open("out_si.dat"); output_ph.open("out_ph.dat"); - output_poro.open("out_poro.ca"); + output_poro.open("out_poro.dat"); + output_saca.open("out_sa_ca.dat"); + output_sasi.open("out_sa_si.dat"); std::vector mineral_out(thedatabase->nb_mineral); for (int mineral: thedatabase->range_mineral()) { mineral_out[mineral].open("out_mineral_"+std::to_string(mineral)+".dat"); } SIASaturatedReactiveTransportSolver solver(themesh, thedatabase, parameters); solver.apply_boundary_conditions(bcs); solver.get_variables().scale_volume(); //solver.use_snia(); solver.use_sia(100); solver.get_options().residual_tolerance = 1e-3; solver.get_options().threshold_update_disable_speciation = 1e-16; solver.get_options().threshold_update_transport = 1e-10; + solver.get_options().good_enough_transport_residual_tolerance = 1e-2; // Transport - solver.get_options().transport_options.residuals_tolerance = 1e-4; - solver.get_options().transport_options.threshold_stationary_point = 1.0; + solver.get_options().transport_options.residuals_tolerance = 1e-12; + //solver.get_options().transport_options.step_tolerance = 1e-10; + solver.get_options().transport_options.threshold_stationary_point = 1e-5; solver.get_options().transport_options.quasi_newton = 5; solver.get_options().transport_options.sparse_solver = SparseSolver::GMRES; + //solver.get_options().transport_options.alpha = 0.0; solver.get_options().transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang; // Speciation + solver.get_options().speciation_options.solver_options.fvectol = 1e-8; solver.get_options().speciation_options.solver_options.use_crashing = false; solver.get_options().speciation_options.solver_options.use_scaling = false; solver.get_options().speciation_options.solver_options.max_iter = 100; solver.get_options().speciation_options.solver_options.maxstep = 50; solver.get_options().speciation_options.solver_options.non_monotone_linesearch = true; + solver.get_options().speciation_options.solver_options.threshold_cycling_linesearch = 1e-8; std::cout << "id SiO2,am : " << dbhandler.mineral_label_to_id("SiO2,am") << std::endl; //thedatabase->stability_mineral( // dbhandler.mineral_label_to_id("SiO2,am")) = database::MineralStabilityClass::cannot_dissolve; //double dt = 5e4; double total=0 ; uindex_t k = 1; - + scalar_t dt = 25.0; while (total < 100*24*3600) { - double dt = get_dt(parameters, themesh); total += dt; + start_again: std::cout << " iterations : " << k << " - dt : " << dt << std::endl; SIASaturatedReactiveTransportSolverReturnCode retcode = solver.solve_timestep(dt); std::cout << " Return code : " << (int) retcode << " - Nb iter : " << solver.get_perfs().nb_iterations << std::endl << " - Nb iter transport : " << solver.get_perfs().nb_iterations_transport << std::endl << " - Nb iter chemistry : " << solver.get_perfs().nb_iterations_chemistry << std::endl; if (retcode != SIASaturatedReactiveTransportSolverReturnCode::Success) + { + std::cout << "Current residual : " << solver.get_perfs().current_residual << std::endl; throw std::runtime_error("Failed to converge : "+std::to_string((int) retcode)); + //dt = dt/1.5; + //goto start_again; + } + + //if (solver.get_perfs().nb_iterations <= 10) dt*=1.1; SIASaturatedVariables& variables = solver.get_variables(); //std::cout << variables.total_mobile_concentrations() << std::endl; // compute porosity for (index_t node: themesh->range_nodes()) { double volume = variables.volume_minerals(node); double cell_volume = themesh->get_volume_cell(node); specmicp_assert(volume < cell_volume); double porosity = (cell_volume - volume*1e6) / cell_volume; parameters->porosity(node) = porosity; parameters->diffusion_coefficient(node) = - porosity>0.92?1e4*std::exp(9.95*0.92-29.08):1e4*std::exp(9.95*porosity-29.08); - //std::cout << parameters->diffusion_coefficient(node); + porosity>0.92?2.219e-5:1e4*std::exp(9.95*porosity-29.08); } - parameters->diffusion_coefficient(0) =1e4*std::exp(9.95*0.92-29.08); // parameters->diffusion_coefficient(1); + std::cout << variables.component_concentration(1, id_ca) << std::endl; + parameters->diffusion_coefficient(0) = 2.219e-5; // parameters->diffusion_coefficient(1); //std::cout << variables.speciation_variables() << std::endl; if (k==1 or (k < 1000 and k% 100 == 0) or k % 1000 == 0) { output_ca << total << " " << total/(3600.0*24.0); output_aqca << total << " " << total/(3600.0*24.0); output_si << total << " " << total/(3600.0*24.0); output_poro << total << " " << total/(3600.0*24.0); output_ph << total << " " << total/(3600.0*24.0); + output_saca << total << " " << total/(3600.0*24.0); + output_sasi << total << " " << total/(3600.0*24.0); for(index_t node: themesh->range_nodes()) { output_ca << " " << variables.nodal_component_update_total_amount(node, id_ca); output_aqca << " " << variables.total_mobile_concentration(node, id_ca); output_si << " " << variables.nodal_component_update_total_amount(node, id_si); output_ph << " " << 14 + variables.log_gamma_component(node, id_oh) + variables.log_component_concentration(node, id_oh); double volume = variables.volume_minerals(node); double cell_volume = themesh->get_volume_cell(node); output_poro << " " << (cell_volume - volume*1e6) / cell_volume; + output_saca << " " << variables.component_concentration_in_minerals(node, id_ca); + output_sasi << " " << variables.component_concentration_in_minerals(node, id_si); } output_ca << std::endl; output_aqca << std::endl; output_si << std::endl; output_poro << std::endl; output_ph << std::endl; + output_saca << std::endl; + output_sasi << std::endl; for (index_t mineral: thedatabase->range_mineral()) { mineral_out[mineral] << total << " " << total/(3600.0*24.0); for(index_t node: themesh->range_nodes()) { mineral_out[mineral] << " " << variables.mineral_amount(node, mineral); } mineral_out[mineral] << std::endl; } } ++k; } output_ca.close(); output_aqca.close(); output_si.close(); output_poro.close(); + output_saca.close(); + output_sasi.close(); for (index_t mineral: thedatabase->range_mineral()) { mineral_out[mineral].close(); } return EXIT_SUCCESS; } int leaching_neutrality() { std::shared_ptr thedatabase = set_database(); Vector compo_oxyde(4), compo_species(4); compo_oxyde << 62.9, 20.6, 5.8, 3.1; compo_from_oxyde(compo_oxyde, compo_species); scalar_t radius = 3.5; //cm scalar_t height = 8; //cm index_t nb_nodes = 36; database::Database dbhandler(thedatabase); //mesh::Mesh1DPtr themesh = mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, height); mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 5); //auto parameters = std::make_shared(nb_nodes, 1e-8, 0.2); std::shared_ptr parameters = std::make_shared( nb_nodes, // nb_nodes, thedatabase->nb_component, // nb_component, thedatabase->nb_aqueous, // nb_aqueous, 2e-5, // the_diffusion_coefficient, 0.25, // porosity 0.5e-3 //reduction_factor ); //parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("H[+]")) = 9e-5; //parameters->diff_coeff_component(dbhandler.component_label_to_id("HO[-]")) = 9e-5; parameters->diff_coeff_secondary(dbhandler.aqueous_label_to_id("Si(OH)4")) = 1e-7; parameters->diff_coeff_component(dbhandler.component_label_to_id("SiO(OH)3[-]")) = 1e-7; parameters->porosity(0) = 1.0; parameters->reduction_factor(0) = 1.0; EquilibriumState initial_state = get_initial_state(thedatabase, compo_species, 0.5); SIABoundaryConditions bcs(nb_nodes); bcs.list_initial_states.push_back(initial_state); bcs.list_initial_states.push_back(get_blank_state(thedatabase)); bcs.bs_types[0] = BoundaryConditionType::FixedComposition; bcs.initial_states[0] = 1; index_t id_ca = dbhandler.component_label_to_id("Ca[2+]"); index_t id_si = dbhandler.component_label_to_id("SiO(OH)3[-]"); std::ofstream output_ca, output_aqca, output_si, output_poro; output_ca.open("out_ca.dat"); output_aqca.open("out_aqca.dat"); output_si.open("out_si.dat"); output_poro.open("out_poro.ca"); std::vector mineral_out(thedatabase->nb_mineral); for (int mineral: thedatabase->range_mineral()) { mineral_out[mineral].open("out_mineral_"+std::to_string(mineral)+".dat"); } SIASaturatedReactiveTransportNeutralitySolver solver(themesh, thedatabase, parameters); solver.apply_boundary_conditions(bcs); solver.get_variables().scale_volume(); solver.use_sia(100); solver.get_options().residual_tolerance = 1e-3; solver.get_options().threshold_update_disable_speciation = 1e-16; solver.get_options().threshold_update_transport = 1e-6; // Transport solver.get_options().transport_options.residuals_tolerance = 1e-5; solver.get_options().transport_options.maximum_iterations = 100; solver.get_options().transport_options.quasi_newton = 2; solver.get_options().transport_options.maximum_step_length = 1000; solver.get_options().transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang; solver.get_options().transport_options.max_iterations_at_max_length = 100; // Speciation solver.get_options().speciation_options.solver_options.use_crashing = false; solver.get_options().speciation_options.solver_options.use_scaling = false; solver.get_options().speciation_options.solver_options.max_iter = 100; solver.get_options().speciation_options.solver_options.maxstep = 50; std::cout << "id SiO2,am : " << dbhandler.mineral_label_to_id("SiO2,am") << std::endl; //thedatabase->stability_mineral( // dbhandler.mineral_label_to_id("SiO2,am")) = database::MineralStabilityClass::cannot_dissolve; //double dt = 5e4; double total=0 ; uindex_t k = 1; std::cout << "Basis :"; for (index_t component: thedatabase->range_component()) { std:: cout << " " << thedatabase->labels_basis[component]; } std::cout << std::endl; std::cout << "Solid :"; for (index_t mineral: thedatabase->range_mineral()) { std:: cout << " " << thedatabase->labels_minerals[mineral]; } std::cout << std::endl; while (total < 10*24*3600) { double dt = 3600/2; total += dt; std::cout << " iterations : " << k << " - dt : " << dt << std::endl; SIASaturatedReactiveTransportSolverReturnCode retcode = solver.solve_timestep(dt); if (retcode != SIASaturatedReactiveTransportSolverReturnCode::Success) throw std::runtime_error("Failed to converge : "+std::to_string((int) retcode)); SIASaturatedVariables& variables = solver.get_variables(); //std::cout << variables.total_mobile_concentrations() << std::endl; // compute porosity for (index_t node: themesh->range_nodes()) { double volume = variables.volume_minerals(node); double cell_volume = themesh->get_volume_cell(node); specmicp_assert(volume < cell_volume); double porosity = (cell_volume - volume*1e6) / cell_volume; parameters->porosity(node) = porosity; parameters->reduction_factor(node) = porosity>0.6?1e4*std::exp(9.95*0.6-29.08)/2e-5:1e4*std::exp(9.95*porosity-29.08)/2e-5; //std::cout << parameters->diffusion_coefficient(node); } parameters->reduction_factor(0) = parameters->reduction_factor(1); //std::cout << variables.speciation_variables() << std::endl; if (k < 2 or k % 100 == 0) { output_ca << total << " " << total/(3600.0*24.0); output_aqca << total << " " << total/(3600.0*24.0); output_si << total << " " << total/(3600.0*24.0); output_poro << total << " " << total/(3600.0*24.0); for(index_t node: themesh->range_nodes()) { output_ca << " " << variables.nodal_component_update_total_amount(node, id_ca); output_aqca << " " << variables.total_mobile_concentration(node, id_ca); output_si << " " << variables.nodal_component_update_total_amount(node, id_si); double volume = variables.volume_minerals(node); double cell_volume = themesh->get_volume_cell(node); output_poro << " " << (cell_volume - volume*1e6) / cell_volume; } output_ca << std::endl; output_aqca << std::endl; output_si << std::endl; output_poro << std::endl; for (index_t mineral: thedatabase->range_mineral()) { mineral_out[mineral] << total << " " << total/(3600.0*24.0); for(index_t node: themesh->range_nodes()) { mineral_out[mineral] << " " << variables.mineral_amount(node, mineral); } mineral_out[mineral] << std::endl; } } ++k; } output_ca.close(); output_aqca.close(); output_si.close(); output_poro.close(); for (index_t mineral: thedatabase->range_mineral()) { mineral_out[mineral].close(); } return EXIT_SUCCESS; } int main() { specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; specmicp::logger::ErrFile::stream() = &std::cerr; return leaching(); //return leaching_neutrality(); }