diff --git a/examples/reactmicp/saturated_react/react_leaching.cpp b/examples/reactmicp/saturated_react/react_leaching.cpp index 1e08314..5f01742 100644 --- a/examples/reactmicp/saturated_react/react_leaching.cpp +++ b/examples/reactmicp/saturated_react/react_leaching.cpp @@ -1,451 +1,449 @@ #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 #include using namespace specmicp; // Initial conditions // ================== specmicp::RawDatabasePtr leaching_database() { specmicp::database::Database thedatabase("../data/cemdata_specmicp.js"); std::map 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::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); 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); + Vector x; AdimensionalSystemSolver solver(raw_data, constraints, options); + solver.initialise_variables(x, 0.75, -4.0); 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; + true_var->upscaling_variables().setZero(); for (index_t node=0; nodenb_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; nodenb_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; nodenb_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 list_initial_composition = {initial, blank}; std::vector index_initial_state(nb_nodes, 0); for (int i=0; i< additional_nodes; ++i) index_initial_state[i] = 1; std::vector 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(constraints, options); + std::make_shared(nb_nodes, constraints, options); std::shared_ptr transport_stagger = std::make_shared(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(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; } } diff --git a/examples/specmicp/adimensional/momas_thermo.cpp b/examples/specmicp/adimensional/momas_thermo.cpp index 0967cce..24ca7d3 100644 --- a/examples/specmicp/adimensional/momas_thermo.cpp +++ b/examples/specmicp/adimensional/momas_thermo.cpp @@ -1,248 +1,249 @@ #include #include "utils/log.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp/adimensional/equilibrium_curve.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional_kinetics/kinetic_variables.hpp" #include "specmicp/adimensional_kinetics/kinetic_model.hpp" #include "specmicp/adimensional_kinetics/kinetic_system_euler_solver.hpp" #include "specmicp/adimensional_kinetics/kinetic_system_solver.hpp" specmicp::RawDatabasePtr get_momas_db() { specmicp::database::Database thedatabase("../data/momas_benchmark.js"); thedatabase.remove_components({thedatabase.component_label_to_id("X5")}); return thedatabase.get_database(); } class AutomaticReactionPath: public specmicp::EquilibriumCurve { public: AutomaticReactionPath() { specmicp::RawDatabasePtr raw_data = get_momas_db(); set_database(raw_data); std::cout << raw_data->nu_aqueous << std::endl; std::cout << raw_data->logk_aqueous << std::endl; std::cout << raw_data->nu_mineral << std::endl; std::cout << raw_data->logk_mineral << std::endl; specmicp::Vector total_concentrations(raw_data->nb_component); total_concentrations << 1, 0.0, -3.0, 0.0, 1.0; constraints() = specmicp::AdimensionalSystemConstraints(total_concentrations); constraints().disable_conservation_water(); constraints().surface_model.model_type = specmicp::SurfaceEquationType::Equilibrium; constraints().surface_model.concentration = 1.0; specmicp::AdimensionalSystemSolverOptions& options = solver_options(); options.solver_options.maxstep = 1.0; options.solver_options.max_iter = 200; options.solver_options.maxiter_maxstep = 200; options.solver_options.use_crashing = false; options.solver_options.use_scaling = false; options.solver_options.non_monotone_linesearch = 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-8; options.solver_options.steptol = 1e-10; options.system_options.non_ideality_tolerance = 1e-10; options.system_options.non_ideality = false; //solve_first_problem(); specmicp::Vector variables; specmicp::AdimensionalSystemSolver solver(raw_data, constraints(), solver_options()); solver.initialise_variables(variables, 0.8, - {{"X2", -3}, {"X4", -10}}, + {{"X2", -3}, {"X4", -11}}, {}, 0.5 ); + solver.run_pcfm(variables); specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, false); std::cout << (int) current_perf.return_code << std::endl; solution_vector() = variables; initialize_solution(solver.get_raw_solution(variables)); std::cout << variables << std::endl; std::cout << "Buffering \t"; for (specmicp::index_t mineral: raw_data->range_mineral()) { std::cout << "\t" << raw_data->labels_minerals[mineral]; } std::cout << "\t S"; for (specmicp::index_t sorbed: raw_data->range_sorbed()) { std::cout << "\t" << raw_data->labels_sorbed[sorbed]; } std::cout << std::endl; output(); } void output() { specmicp::AdimensionalSystemSolutionExtractor sol(current_solution(), database(), solver_options().units_set); std::cout << constraints().total_concentrations(1); for (specmicp::index_t mineral: database()->range_mineral()) { std::cout << "\t" << sol.mole_concentration_mineral(mineral); } std::cout << "\t" << sol.free_surface_concentration(); for (specmicp::index_t sorbed: database()->range_sorbed()) { std::cout << "\t" << sol.sorbed_species_molalities(sorbed); } std::cout << std::endl; } void update_problem() { constraints().total_concentrations(1) += 0.1; constraints().total_concentrations(2) += 0.1; constraints().total_concentrations(3) += 0.1; } private: specmicp::index_t id_h2o; specmicp::index_t id_ho; specmicp::index_t id_hco3; specmicp::index_t id_co2g; specmicp::index_t id_ca; }; class MomasKinetics: public specmicp::kinetics::AdimKineticModel { public: MomasKinetics(): specmicp::kinetics::AdimKineticModel({0}), m_id_cc(0), m_id_c3(2), m_id_x4(4) {} //! \brief Compute the kinetic rates and store them in dydt void compute_rate( specmicp::scalar_t t, const specmicp::Vector& y, specmicp::kinetics::AdimKineticVariables& variables, specmicp::Vector& dydt ) { dydt.resizeLike(y); specmicp::AdimensionalSystemSolution& solution = variables.equilibrium_solution(); specmicp::scalar_t factor = std::pow(1e3*0.5*solution.secondary_molalities(m_id_c3),3)/ std::pow(1e3*0.5*pow10(solution.main_variables(m_id_x4)), 2); factor = 0.2*factor-1.0; specmicp::scalar_t k_constant = 10.0; if (factor >= 0) k_constant = 1e-2; dydt(0) = k_constant*factor; } private: specmicp::index_t m_id_cc; specmicp::index_t m_id_c3; specmicp::index_t m_id_x4; }; void solve_kinetics_problem() { specmicp::RawDatabasePtr raw_data = get_momas_db(); specmicp::Vector total_concentrations(6); total_concentrations << 1, 0.6, -2.4, 0.6, 1.0, 1.0; specmicp::AdimensionalSystemConstraints constraints(total_concentrations); constraints.disable_conservation_water(); 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.non_monotone_linesearch = false; options.solver_options.factor_descent_condition = -1; options.solver_options.factor_gradient_search_direction = 50; options.solver_options.projection_min_variable = 1e-9; options.solver_options.fvectol = 1e-8; options.solver_options.steptol = 1e-10; options.system_options.non_ideality_tolerance = 1e-10; options.system_options.non_ideality = false; specmicp::Vector equil_variables(raw_data->nb_component+raw_data->nb_mineral); equil_variables << 0.8, -5, -1, -5, -1, -3, 0.0; specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(equil_variables, false); std::cout << (int) current_perf.return_code << std::endl; specmicp::AdimensionalSystemSolution equil_solution = solver.get_raw_solution(equil_variables); std::shared_ptr model = std::make_shared(); specmicp::Vector mineral_concentrations(1); mineral_concentrations << 5.0; specmicp::kinetics::AdimKineticSystemEulerSolver kinetic_solver( //specmicp::kinetics::AdimKineticSystemSolver kinetic_solver( model, total_concentrations, mineral_concentrations, constraints, equil_solution, raw_data); options.solver_options.use_scaling = false; kinetic_solver.get_options().speciation_options = options; kinetic_solver.solve(0.005, 1.0); std::cout << "Final amount : " << kinetic_solver.variables().concentration_mineral(0) << std::endl; } int main() { specmicp::logger::ErrFile::stream() = &std::cerr; - specmicp::stdlog::ReportLevel() = specmicp::logger::Warning; + specmicp::stdlog::ReportLevel() = specmicp::logger::Debug; // solve_kinetics_problem(); AutomaticReactionPath test_automatic; for (int i=0; i<20; ++i) { test_automatic.run_step(); } return EXIT_SUCCESS; }