diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index a35c686..d306622 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,72 +1,79 @@ set(PROJECT_EXAMPLE_DIR ${CMAKE_CURRENT_SOURCE_DIR}) # ========== SpecMiCP ===================================================== # ---------- adimensional system ------------------------------------------ # thermocarbo add_executable(ex_adim_thermocarbo ${PROJECT_EXAMPLE_DIR}/specmicp/adimensional/thermocarbo.cpp ) target_link_libraries(ex_adim_thermocarbo specmicp specmicp_database) # Momas thermodynamic example add_executable(ex_adim_momas ${PROJECT_EXAMPLE_DIR}/specmicp/adimensional/momas_thermo.cpp ) target_link_libraries(ex_adim_momas specmicp specmicp_database) # equilibrium curve add_executable(ex_adim_equilibriumcurve ${PROJECT_EXAMPLE_DIR}/specmicp/adimensional/equilibrium_curve.cpp ) target_link_libraries(ex_adim_equilibriumcurve specmicp specmicp_database) # ========== ReactMiCP ==================================================== # ---------- leaching ----------------------------------------------------- add_executable(leaching ${PROJECT_EXAMPLE_DIR}/reactmicp/leaching.cpp ) target_link_libraries(leaching reactmicp specmicp specmicp_database) # ---------- carbonation -------------------------------------------------- add_executable(carbonation ${PROJECT_EXAMPLE_DIR}/reactmicp/carbonation.cpp ) target_link_libraries(carbonation reactmicp specmicp specmicp_database) # ---------- diffusion ---------------------------------------------------- add_executable(diffusion ${PROJECT_EXAMPLE_DIR}/reactmicp/diffusion.cpp ) target_link_libraries(diffusion reactmicp specmicp specmicp_database) # ---------- EquilibriumCurve ---------------------------------------------------- add_executable(ex_reactmicp_equilibriumcurve ${PROJECT_EXAMPLE_DIR}/reactmicp/equilibrium_curve.cpp ) target_link_libraries(ex_reactmicp_equilibriumcurve reactmicp dfpm specmicp specmicp_database) # New flexible reactive transport solver # ======================================= # ---------- Leaching ---------------------------------------------------- add_executable(ex_reactmicp_leaching ${PROJECT_EXAMPLE_DIR}/reactmicp/saturated_react/react_leaching.cpp ) target_link_libraries(ex_reactmicp_leaching reactmicp specmicp specmicp_database) + +# ---------- Momas benchmark --------------------------------------------- + +add_executable(ex_reactmicp_momas + ${PROJECT_EXAMPLE_DIR}/reactmicp/saturated_react/momas_benchmark.cpp +) +target_link_libraries(ex_reactmicp_momas reactmicp specmicp specmicp_database) diff --git a/examples/reactmicp/saturated_react/momas_benchmark.cpp b/examples/reactmicp/saturated_react/momas_benchmark.cpp new file mode 100644 index 0000000..fc3dc8e --- /dev/null +++ b/examples/reactmicp/saturated_react/momas_benchmark.cpp @@ -0,0 +1,525 @@ +#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 "specmicp/adimensional/equilibrium_curve.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 "dfpm/meshes/uniform_mesh1d.hpp" + +#include "utils/log.hpp" +#include +#include + +#include "database/selector.hpp" + +using namespace specmicp; + +RawDatabasePtr get_momas_db() +{ + database::Database thedatabase("../data/momas_benchmark.js"); + return thedatabase.get_database(); +} + +void prepare_db_for_easy_case(RawDatabasePtr raw_data) +{ + database::DatabaseSelector selector(raw_data); + selector.remove_component({selector.component_label_to_id("X5"),}); + specmicp_assert(raw_data->nb_component == 5); + selector.remove_all_minerals(); + //std::cout << selector.aqueous_label_to_id("C6") << " - " << selector.aqueous_label_to_id("C7") << std::endl; + selector.remove_aqueous({selector.aqueous_label_to_id("C6"), selector.aqueous_label_to_id("C7")}); + //specmicp_assert(raw_data->nb_aqueous == 5); +} + + +class AutomaticReactionPath: public specmicp::EquilibriumCurve +{ +public: + + AutomaticReactionPath(specmicp::RawDatabasePtr raw_data) + { + set_database(raw_data); + for (index_t aqueous: raw_data->range_aqueous()) + { + std::cout << raw_data->labels_aqueous[aqueous] << std::endl; + } + + 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; + std::cout << raw_data->nu_sorbed << std::endl; + std::cout << raw_data->logk_sorbed << std::endl; + + specmicp::Vector total_concentrations(raw_data->nb_component); + total_concentrations << 1, 0, -2.0, 0, 2.0; + + constraints() = specmicp::AdimensionalSystemConstraints(0.25*total_concentrations); + constraints().disable_conservation_water(); + constraints().inert_volume_fraction = 0.75; + constraints().surface_model.model_type = specmicp::SurfaceEquationType::Equilibrium; + constraints().surface_model.concentration = 1.0; + + specmicp::AdimensionalSystemSolverOptions& options = solver_options(); + options.solver_options.maxstep = 0.5; + 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 = 100; + options.solver_options.projection_min_variable = 1e-9; + options.solver_options.fvectol = 1e-10; + options.solver_options.steptol = 1e-12; + options.system_options.non_ideality_tolerance = 1e-10; + options.system_options.non_ideality = false; + //options.units_set.length = units::LengthUnit::decimeter; + + //solve_first_problem(); + + specmicp::Vector variables; + specmicp::AdimensionalSystemSolver solver(raw_data, constraints(), solver_options()); + + solver.initialise_variables(variables, + 0.25, + {{"X2", -2}, {"X4", -6}}, + {}, + -0.3 + ); + solver.run_pcfm(variables); + std::cout << variables << std::endl; + + specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, false); + std::cout << (int) current_perf.return_code << std::endl; + + specmicp_assert((int) current_perf.return_code == 2); + + 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() + { + //std::cout << "\n --------- \n New problem \n --------- \n"; + //constraints().total_concentrations(1) += 0.1; + constraints().total_concentrations(2) += -0.15; + //constraints().total_concentrations(3) += 0.1; + constraints().total_concentrations(4) += 0.15; + constraints().surface_model.concentration += 0; + + } + +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; +}; + +AdimensionalSystemSolution easy_material_a( + RawDatabasePtr raw_data, + const AdimensionalSystemConstraints& constraints, + const AdimensionalSystemSolution& initial_solution, + const AdimensionalSystemSolverOptions& options + ) +{ + specmicp::Vector variables = initial_solution.main_variables; + std::cout << constraints.total_concentrations << std::endl; + specmicp::AdimensionalSystemSolver solver(raw_data, constraints, initial_solution, options); + specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, false); + std::cout << (int) current_perf.return_code << std::endl; + + AdimensionalSystemSolutionExtractor sol(solver.get_raw_solution(variables), + raw_data, + options.units_set); + + std::cout << sol.total_aqueous_concentration(2) << " - " << sol.total_immobile_concentration(2) << std::endl; + std::cout << sol.density_water() << std::endl; + std::cout << 0.25*sol.total_aqueous_concentration(2) + sol.total_immobile_concentration(2) + << " =? " + << -2*0.25 << std::endl; + return solver.get_raw_solution(variables); +} + + +AdimensionalSystemSolution easy_material_b( + RawDatabasePtr raw_data, + const AdimensionalSystemConstraints& constraints, + const AdimensionalSystemSolution& initial_solution, + const specmicp::AdimensionalSystemSolverOptions& options + ) +{ + specmicp::Vector variables = initial_solution.main_variables; + + specmicp::AdimensionalSystemSolver solver(raw_data, constraints, initial_solution, options); + specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, false); + std::cout << (int) current_perf.return_code << std::endl; + + return solver.get_raw_solution(variables); +} + +AdimensionalSystemSolution easy_bc( + RawDatabasePtr raw_data, + const AdimensionalSystemConstraints& constraints, + const specmicp::AdimensionalSystemSolverOptions& options + ) +{ + specmicp::Vector variables; + specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options); + solver.initialise_variables(variables, 1.0, -4.0); + solver.run_pcfm(variables); + specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, false); + std::cout << (int) current_perf.return_code << std::endl; + AdimensionalSystemSolutionExtractor ext(solver.get_raw_solution(variables), + raw_data, + options.units_set); + scalar_t aq_concentration = ext.total_aqueous_concentration(2); + std::cout << "Aqueous concentration 2 : " << aq_concentration << std::endl; + return solver.get_raw_solution(variables); +} + +using namespace specmicp::reactmicp; +using namespace specmicp::reactmicp::solver; +using namespace specmicp::reactmicp::systems::satdiff; + +class MoMasUspcaling: public solver::UpscalingStaggerBase +{ +public: + + MoMasUspcaling(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); + + true_var->upscaling_variables().setZero(); + for (index_t node=0; nodenb_nodes(); ++node) + { + // if material b + const scalar_t coord = true_var->get_mesh()->get_position(node); + if (coord < 1.1 and coord >= 1.0) + { + true_var->porosity(node) = 0.5; + true_var->permeability(node) = 1e-5; + true_var->diffusion_coefficient(node) = 3.3e-4; + true_var->fluid_velocity(node) = 5.5e-3; + } + else // material a + { + true_var->porosity(node) = 0.25; + true_var->permeability(node) = 1e-2; + true_var->diffusion_coefficient(node) = 5.5e-5; + true_var->fluid_velocity(node) = 5.5e-3; + } + std::cout << "node : " << node << " - porosity : " << true_var->porosity(node) << std::endl; + } + } + + //! \brief Initialize the stagger at the beginning of an iteration + void initialize_timestep(scalar_t dt, VariablesBasePtr var) + { + } + + //! \brief Solve the equation for the timestep + StaggerReturnCode restart_timestep(VariablesBasePtr var) + { + + return StaggerReturnCode::ResidualMinimized; + } + +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 database = get_momas_db(); + prepare_db_for_easy_case(database); + + units::UnitsSet the_units = units::UnitsSet(); + the_units.length = units::LengthUnit::decimeter; + + + specmicp::Vector total_concentrations(database->nb_component); + total_concentrations << 1, 0.0, -2.0, 0.0, 2.0; + specmicp::AdimensionalSystemConstraints constraints_a(0.25*total_concentrations); + constraints_a.disable_conservation_water(); + constraints_a.surface_model.model_type = SurfaceEquationType::Equilibrium; + constraints_a.surface_model.concentration = 1.0; + constraints_a.inert_volume_fraction = 0.75; + + specmicp::AdimensionalSystemConstraints constraints_b(0.5*total_concentrations); + constraints_b.disable_conservation_water(); + constraints_b.surface_model.model_type = SurfaceEquationType::Equilibrium; + constraints_b.surface_model.concentration = 10.0; + constraints_b.inert_volume_fraction = 0.5; + + specmicp::Vector total_concentrations_bc(database->nb_component); + total_concentrations_bc << 1, 0.3, 0.3, 0.3, 0.0; + specmicp::AdimensionalSystemConstraints constraints_bc(total_concentrations_bc); + constraints_bc.disable_conservation_water(); + constraints_bc.surface_model.model_type = SurfaceEquationType::NoEquation; + constraints_bc.surface_model.concentration = 0.0; + constraints_bc.inert_volume_fraction = 0.0; + + + AdimensionalSystemSolverOptions 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-6; + options.solver_options.steptol = 1e-14; + options.solver_options.threshold_cycling_linesearch = 1e-8; + options.system_options.non_ideality_tolerance = 1e-10; + options.system_options.non_ideality = false; + options.system_options.restart_concentration = -2; + options.use_pcfm = true; + options.units_set = the_units; + // options.force_pcfm = true; + + AdimensionalSystemSolution mat_bc = easy_bc(database, constraints_bc, options); + + + + AutomaticReactionPath test_automatic(database); + + for (int i=0; i<10; ++i) + { + test_automatic.run_step(); + } + AdimensionalSystemSolution pre_solution = test_automatic.current_solution(); + AdimensionalSystemSolution mat_a = easy_material_a(database, constraints_a, pre_solution, options); + AdimensionalSystemSolution mat_b = easy_material_b(database, constraints_b, pre_solution, options); + + options.solver_options.maxstep = 1.0; + options.solver_options.use_scaling = true; + options.solver_options.non_monotone_linesearch = true; + //options.force_pcfm = true; + + scalar_t tot_length = 2.1; + scalar_t dx = 0.01; + + index_t nb_nodes = tot_length/dx +1; + + mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d(nb_nodes, dx, 1.0); + + std::vector list_constraints = {constraints_a, + constraints_b, + constraints_bc}; + std::vector list_initial_composition = {mat_a, mat_b, mat_bc}; + std::vector index_initial_state(nb_nodes, 0); + std::vector index_constraints(nb_nodes, 0); + + index_t start_mat_b = 1.0/dx; + index_t end_mat_b = 1.1/dx; + + for (index_t node=start_mat_b; node list_fixed_nodes = {0, }; + std::map list_slaves_nodes = {{nb_nodes-1, nb_nodes-2}}; + std::vector list_immobile_components = {0}; + + systems::satdiff::SaturatedVariablesPtr variables = + systems::satdiff::init_variables(the_mesh, database, the_units, + list_fixed_nodes, + list_initial_composition, + index_initial_state); + variables->displacement()(1) = 0.3; + variables->displacement()(2) = 0.3; + variables->displacement()(3) = 0.3; + variables->displacement()(6) = 0.0; + variables->displacement()(7) = 0.0; + variables->displacement()(8) = 0.0; + + std::cout << variables->displacement() << std::endl; + + ChemistryStaggerPtr equilibrium_stagger = + std::make_shared(list_constraints, index_constraints, options); + + std::shared_ptr transport_stagger = + std::make_shared( + variables, list_fixed_nodes, list_slaves_nodes, list_immobile_components); + + 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-12; + 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(database, 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 = 50; //500; + solver.get_options().residuals_tolerance = 1e-2; + + solver.get_options().absolute_residuals_tolerance = 1e-16; + solver.get_options().implicit_upscaling = false; + + +// std::cout << variables->displacement() << std::endl; +// std::cout << "Final results " << std::endl; +// transport_stagger->initialize_timestep(0.1, variables); +// transport_stagger->restart_timestep(variables); + +// std::cout << variables->displacement() << std::endl; + + //std::cout << variables->upscaling_variables() << std::endl; + scalar_t total = 0.0; + scalar_t dt=0.001; + + + //std::cout << variables->displacement() << std::endl; + + std::ofstream output_x1, output_x2, output_x3, output_x4; + std::ofstream output_sx1, output_sx2, output_sx3, output_sx4; + std::ofstream output_c1, output_s; + output_x1.open("out_momas_x1.dat"); + output_x2.open("out_momas_x2.dat"); + output_x3.open("out_momas_x3.dat"); + output_x4.open("out_momas_x4.dat"); + output_sx1.open("out_momas_sx1.dat"); + output_sx2.open("out_momas_sx2.dat"); + output_sx3.open("out_momas_sx3.dat"); + output_sx4.open("out_momas_sx4.dat"); + output_c1.open("out_momas_c1.dat"); + output_s.open("out_momas_s.dat"); + index_t cnt = 0; + + while (total < 20000*dt) + { + solver.solve_timestep(dt, variables); + //std::cout<< variables->displacement() << std::endl; + total += dt; + + std::cout << "Time : " << total << std::endl + << "Iter : " << solver.get_perfs().nb_iterations << std::endl + << "Residuals : " << solver.get_perfs().residuals << std::endl; + + if (cnt % 10) + { + output_x1 << total; + output_x2 << total; + output_x3 << total; + output_x4 << total; + output_sx1 << total; + output_sx2 << total; + output_sx3 << total; + output_sx4 << total; + output_c1 << total; + output_s << total; + for (index_t node: the_mesh->range_nodes()) + { + output_x1 << "\t" << variables->aqueous_concentration(node, 1, variables->displacement()); + output_x2 << "\t" << variables->aqueous_concentration(node, 2, variables->displacement()); + output_x3 << "\t" << variables->aqueous_concentration(node, 3, variables->displacement()); + output_x4 << "\t" << variables->aqueous_concentration(node, 4, variables->displacement()); + output_sx1 << "\t" << variables->solid_concentration(node, 1, variables->displacement()); + output_sx2 << "\t" << variables->solid_concentration(node, 2, variables->displacement()); + output_sx3 << "\t" << variables->solid_concentration(node, 3, variables->displacement()); + output_sx4 << "\t" << variables->solid_concentration(node, 4, variables->displacement()); + output_c1 << "\t" << variables->equilibrium_solution(node).secondary_molalities(0); + output_s << "\t" << pow10(variables->equilibrium_solution(node).main_variables(5)); + } + output_x1 << std::endl; + output_x2 << std::endl; + output_x3 << std::endl; + output_x4 << std::endl; + output_sx1 << std::endl; + output_sx2 << std::endl; + output_sx3 << std::endl; + output_sx4 << std::endl; + output_c1 << std::endl; + output_s << std::endl; + } + ++ cnt; + } + return EXIT_SUCCESS; +}