diff --git a/examples/reactmicp/saturated_react/carbonation.cpp b/examples/reactmicp/saturated_react/carbonation.cpp
index 09333e1..a231396 100644
--- a/examples/reactmicp/saturated_react/carbonation.cpp
+++ b/examples/reactmicp/saturated_react/carbonation.cpp
@@ -1,510 +1,509 @@
 // ======================================================
 
 //    Leaching of a cement paste in CO2-saturated water
 
 // =======================================================
 
 // This file can be used to learn how to use ReactMiCP
 
 
 // Include ReactMiCP
 #include "reactmicp.hpp"
 
 // Other includes that would be needed later :
 #include <Eigen/QR>
 
 #include "physics/laws.hpp"
 #include "utils/timer.hpp"
 #include <functional>
 #include <iostream>
 #include <fstream>
 
 
 // we use the following namespaces
 using namespace specmicp;
 using namespace reactmicp;
 using namespace reactmicp::systems;
 
 using namespace specmicp::reactmicp::systems::satdiff;
 using namespace reactmicp::solver;
 
 // We declare some functions that we will be using :
 
 // Boundary and initial conditions
 // ===============================
 
 // Return the chemistry constraints for the inital condition
 AdimensionalSystemConstraints get_cement_initial_constraints(
         RawDatabasePtr the_database,
         const units::UnitsSet& units_set);
 // return the chemistry constraints for the boundary condition
 AdimensionalSystemConstraints get_bc_initial_constraints(
         RawDatabasePtr the_database,
         const units::UnitsSet& units_set);
 
 
 Vector initial_compo(const Vector& publi);
 
 void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species);
 
 // The mesh
 // ========
 
 // There is two version
 // a uniform mesh (cf config)
 // and a non uniform mesh
 mesh::Mesh1DPtr get_mesh();
 
 // ===================
 //       MAIN
 // ===================
 
 int main()
 {
     // initialize eigen
     // This is needed if OpenMP is used
     Eigen::initParallel();
 
     // Setup the log
     // The output will be on cerr
     specmicp::logger::ErrFile::stream() = &std::cerr;
     // And only Warning and error messages will be printed
     specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
 
 
     // Configuration
 
     auto config = specmicp::io::parse_config_file("carbonation.yaml");
 
     // The database
     // -------------
 
     // Obtain the database
     RawDatabasePtr the_database = specmicp::io::configure_database(config);
 
     // SpecMiCP options
     // ----------------
 
     // Set the options
     AdimensionalSystemSolverOptions specmicp_options;
     specmicp::io::configure_specmicp_options(specmicp_options, config);
     units::UnitsSet units_set =  specmicp_options.units_set;
 
 
     // The mesh
     // --------
 
     // Obtain the mesh
     mesh::Mesh1DPtr the_mesh = specmicp::io::configure_mesh(config);
     index_t nb_nodes = the_mesh->nb_nodes(); // just boilerplate to make things easier
 
     // Print the mesh
     io::print_mesh("out_carbo_mesh.dat", the_mesh, units::LengthUnit::centimeter);
 
 
     // Initial and boundary conditions
     // --------------------------------
 
     // Obtain the initial constraints
     // Note : the database is modified at this point
     // Only the relevant components are kept at this point
     AdimensionalSystemConstraints cement_constraints = get_cement_initial_constraints(the_database, units_set);
 
     // Freeze the database
     the_database->freeze_db();
 
 
     // Obtain the initial condition
     AdimensionalSystemSolution cement_solution = solve_equilibrium(
                 the_database,
                 cement_constraints,
                 specmicp_options);
 
     if(not the_database->is_valid())
         throw std::runtime_error("The database has been changed during the initial condition computation");
 
 
     // Just to check that everything is ok !
     AdimensionalSystemSolutionExtractor extractor(cement_solution, the_database, units_set);
     std::cout << "Porosity : " << extractor.porosity() << std::endl;
     std::cout << "Water volume fraction " << extractor.volume_fraction_water() << std::endl;
 
 
 
     //return EXIT_SUCCESS;
 
     // Obtain the boundary conditions
     AdimensionalSystemConstraints bc_constraints = get_bc_initial_constraints(
                 the_database,
                 specmicp_options.units_set
                 );
     AdimensionalSystemSolution bc_solution = solve_equilibrium(
                 the_database,
                 bc_constraints,
                 specmicp_options
                 );
 
     if(not the_database->is_valid())
         throw std::runtime_error("The database has been changed during the boundary conditions computation");
 
 
     // The constraints
     // --------------
 
     // Boundary condition
     // 'list_fixed_nodes' lists the node with fixed composition
     std::vector<index_t> list_fixed_nodes = {0};
 
 
     // list_constraints lists the different speciation constraints in the system
     // Note : the total concentrations are computed from the initial solution,
     //  and thus different total concentrations do not constitute a reason to create different constraints
     std::vector<specmicp::AdimensionalSystemConstraints> list_constraints = {cement_constraints,
                                                                              bc_constraints};
 
     // index_constraints links a node to the set of constraints that will be used
     std::vector<int> index_constraints(nb_nodes, 0);
     index_constraints[0] = 1; // Boundary conditions
 
     // This is the list of initial composition
     std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition = {cement_solution, bc_solution};
     // This list links a node to its initial conditions
     std::vector<int> index_initial_state(nb_nodes, 0);
     index_initial_state[0] = 1; // boundary condition
 
     // Variables
     // -----------
 
     // Create the main set of variables
     // They will be initialized using the list of initial composition
     satdiff::SaturatedVariablesPtr variables =
             satdiff::init_variables(the_mesh, the_database, units_set,
                                                         list_fixed_nodes,
                                                         list_initial_composition,
                                                         index_initial_state);
 
     // Staggers
     // ---------
 
     // This section initialize the staggers
 
     // The equilibrium stagger
     // - - - - - - - - - - - -
 
     std::shared_ptr<EquilibriumStagger> equilibrium_stagger =
             std::make_shared<satdiff::EquilibriumStagger>(list_constraints, index_constraints, specmicp_options);
 
     // The transport stagger
     // - - - - - - - - - - -
 
     std::vector<index_t> list_fixed_component {0, 1, the_database->get_id_component("Si(OH)4")};
     std::map<index_t, index_t> slave_nodes {};
 
 
     std::shared_ptr<SaturatedTransportStagger> transport_stagger =
             std::make_shared<satdiff::SaturatedTransportStagger>(
                 variables, list_fixed_nodes, slave_nodes, list_fixed_component
                 );
 
     specmicp::io::configure_transport_options(transport_stagger->options_solver(), config);
 
     // The upscaling stagger
     // - - - - - - - - - - -
 
     CappedPowerLawParameters upscaling_param;
     upscaling_param.d_eff_0 = 2.726e-12*1e4;
     upscaling_param.cap = 1e10*1e6;
     upscaling_param.exponent = 3.32;
     upscaling_param.porosity_0 = 0.25;
     upscaling_param.porosity_res = 0.02;
 
     CappedPowerLaw upscaling_law(upscaling_param);
 
     UpscalingStaggerPtr upscaling_stagger =
             std::make_shared<DiffusiveUpscalingStagger>(upscaling_law.get_law(),the_database, units_set);
 
     // Solver
     // ------
 
     // This section initialize the reactive transport solver
 
     ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger);
     specmicp::io::configure_reactmicp_options(solver.get_options(), config);
     transport_stagger->initialize(variables);
     equilibrium_stagger->initialize(variables);
     upscaling_stagger->initialize(variables);
 
 
     // Some basic information about the simulation
     SimulationInformation simul_info = specmicp::io::configure_simulation_information(config);
 
     // Runner
     // -------
 
     // Create the runner : loop through the timestep
 
     specmicp::io::check_mandatory_yaml_node(config, "run", "__main__");
     scalar_t min_dt = specmicp::io::get_yaml_mandatory<scalar_t>(config["run"], "minimum_dt", "run");
     scalar_t max_dt = specmicp::io::get_yaml_mandatory<scalar_t>(config["run"], "maximum_dt", "run");
     scalar_t total = specmicp::io::get_yaml_mandatory<scalar_t>(config["run"], "total_execution_time", "run");
     ReactiveTransportRunner runner(solver, min_dt, max_dt, simul_info);
 
     // Output
     // ------
 
     io::OutputNodalVariables output_policy(the_database, the_mesh, specmicp_options.units_set);
     io::configure_reactmicp_output(output_policy, simul_info, config);
 
     // and bind it to the runner
     runner.set_output_policy(output_policy.get_output_for_reactmicp());
 
     // Solve the problem
     // -----------------
 
     runner.run_until(total, cast_var_to_base(variables));
 
 
     if(not the_database->is_valid())
         throw std::runtime_error("The database has been changed during the computation");
 
     AdimensionalSystemSolutionSaver saver(the_database);
     saver.save_solutions(simul_info.complete_filepath("solutions", "dat"), variables->equilibrium_solutions());
 
     // output information at the end of the timestep
 
     io::print_minerals_profile(the_database, variables, the_mesh, specmicp_options.units_set,
                                simul_info.complete_filepath("profiles_end", "dat"));
     io::print_components_total_aqueous_concentration(the_database, variables,
                                                      the_mesh, specmicp_options.units_set,
                                                      simul_info.complete_filepath("c_profiles_end", "dat"));
     io::print_components_total_solid_concentration(the_database, variables,
                                                    the_mesh, specmicp_options.units_set,
                                                    simul_info.complete_filepath("s_profiles_end", "dat"));
 
     return EXIT_SUCCESS;
 }
 
 AdimensionalSystemConstraints get_cement_initial_constraints(
         RawDatabasePtr the_database,
         const units::UnitsSet& units_set)
 {
 
     Vector oxyde_compo(4);
     oxyde_compo << 66.98, 21.66, 7.19, 2.96;
     Vector species_compo;
     compo_from_oxyde(oxyde_compo, species_compo);
 
     scalar_t mult = 1e-6;
     //species_compo *= 0.947e-2; // poro 0.47
     //species_compo *= 1.08e-2; // poro 0.4
     //species_compo *= 1.25e-2;  // poro 0.3
     species_compo *= 1.1e-2;
     Formulation formulation;
     formulation.mass_solution = 1.0;
     formulation.amount_minerals = {
         {"C3S", species_compo(0)},
         {"C2S", species_compo(1)},
         {"C3A", species_compo(2)},
         {"Gypsum", species_compo(3)},
         {"Calcite", species_compo(0)*1e-3}
 
     };
     formulation.extra_components_to_keep = {"HCO3[-]"};
     formulation.concentration_aqueous =
     {
-        {"Na[+]", mult*1.0298},
-        {"K[+]",  mult*0.0801},
-        {"Cl[-]", mult*0.0001},
-        {"HO[-]", mult*1.8298}
+        {"NaOH", mult*1.0298},
+        {"KOH",  mult*0.0801},
+        {"Cl[-]", mult*0.0001}
     };
     Dissolver dissolver(the_database);
     dissolver.set_units(units_set);
 
     AdimensionalSystemConstraints constraints;
-    constraints.set_charge_keeper(1);
+    constraints.set_charge_keeper(the_database->get_id_component("HO[-]"));
 
     constraints.total_concentrations =  dissolver.dissolve(formulation);
 
     std::cout << constraints.total_concentrations << std::endl;
 
     constraints.set_saturated_system();
 
     return constraints;
 }
 
 AdimensionalSystemConstraints get_bc_initial_constraints(
         RawDatabasePtr the_database,
         const units::UnitsSet& units_set)
 {
     const scalar_t rho_w = laws::density_water(units::celsius(25.0), units_set.length, units_set.mass);
     const scalar_t nacl = 0.3*rho_w;
     const scalar_t kcl = 0.28*rho_w;
 
     Vector total_concentrations;
     total_concentrations.setZero(the_database->nb_component());
     total_concentrations(the_database->get_id_component("K[+]")) = kcl;
     total_concentrations(the_database->get_id_component("Cl[-]")) = nacl + kcl;
     total_concentrations(the_database->get_id_component("Na[+]")) = nacl;
     //total_concentrations(the_database->get_id_component("Si(OH)4")) = 0.0001*rho_w;
 
     AdimensionalSystemConstraints constraints;
     constraints.add_fixed_activity_component(the_database->get_id_component("HO[-]"), -10.3);
     constraints.set_charge_keeper(the_database->get_id_component("HCO3[-]"));
     constraints.total_concentrations =  total_concentrations;
 
     constraints.set_saturated_system();
     constraints.disable_surface_model();
     constraints.fixed_fugacity_cs ={};
 
     return constraints;
 }
 
 AdimensionalSystemSolution get_bc_initial_compo(
         RawDatabasePtr the_database,
         const AdimensionalSystemConstraints& constraints,
         const AdimensionalSystemSolverOptions& options
         )
 {
     Vector variables;
     AdimensionalSystemSolver solver(the_database,
                                     constraints,
                                     options);
 
     micpsolver::MiCPPerformance perf = solver.solve(variables, true);
 
     if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet)
         throw std::runtime_error("Failed to solve initial composition");
 
     std::cout << variables << std::endl;
 
     return solver.get_raw_solution(variables);
 }
 
 void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species)
 {
     constexpr double M_CaO = 56.08;
     constexpr double M_SiO2 = 60.09;
     constexpr double M_Al2O3 = 101.96;
     constexpr double M_SO3 = 80.06;
 
     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::ColPivHouseholderQR<Eigen::MatrixXd> solver(compo_in_oxyde);
     compo_species = solver.solve(compo_oxyde);
 
 }
 
 
 Vector initial_compo(const Vector& publi)
 {
     Matrix publi_stochio(5,5);
     publi_stochio << 1, 0, 0, 0, 0,
                      9, 6, 0, 0, 0,
                      3, 0, 1, 1, 0,
                      3, 0, 1, 3, 0,
                      0, 0, 0, 0, 1;
 
     Matrix cemdata_stochio(5, 5);
     cemdata_stochio << 1, 0, 0, 0, 0,
                        1.0/0.6, 1.0, 0, 0, 0,
                        3, 0, 1, 1, 0,
                        3, 0, 1, 3, 0,
                        0, 0, 0, 0, 1;
 
     Vector oxyde;
     Eigen::ColPivHouseholderQR<Matrix> solver(publi_stochio);
     oxyde = solver.solve(publi);
 
     Vector for_cemdata =  cemdata_stochio*oxyde;
     return for_cemdata;
 }
 
 mesh::Mesh1DPtr get_mesh()
 {
     Vector radius(71);
 
     const scalar_t height = 20;
 
     radius <<
               3.751,
             3.750,
             3.745,
             3.740,
             3.735,
             3.730,
             3.720,
             3.710,
             3.700,
             3.690,
             3.680,
             3.670,
             3.660,
             3.645,
             3.630,
             3.615,
             3.600,
             3.580,
             3.560,
             3.540,
             3.520,
             3.500,
             3.475,
             3.450,
             3.425,
             3.400,
             3.375,
             3.350,
             3.325,
             3.300,
             3.275,
             3.250,
             3.225,
             3.200,
             3.175,
             3.150,
             3.125,
             3.100,
             3.050,
             3.025,
             3.000,
             2.950,
             2.900,
             2.850,
             2.800,
             2.750,
             2.700,
             2.650,
             2.600,
             2.550,
             2.500,
             2.450,
             2.400,
             2.350,
             2.300,
             2.250,
             2.200,
             2.150,
             2.100,
             2.050,
             2.000,
             1.950,
             1.900,
             1.850,
             1.800,
             1.750,
             1.700,
             1.650,
             1.600,
             1.550,
             1.500
             ;
     radius = radius/10;
     return mesh::axisymmetric_mesh1d(radius, height);
 }