diff --git a/tests/reactmicp/systems/saturated_react/equilibrium_stagger.cpp b/tests/reactmicp/systems/saturated_react/equilibrium_stagger.cpp
index 1db8697..a92efca 100644
--- a/tests/reactmicp/systems/saturated_react/equilibrium_stagger.cpp
+++ b/tests/reactmicp/systems/saturated_react/equilibrium_stagger.cpp
@@ -1,124 +1,123 @@
 #include "catch.hpp"
 
 #include "speciation_system.hpp"
 #include "specmicp/adimensional/adimensional_system_solution.hpp"
 #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
 
 #include "specmicp/adimensional/adimensional_system_solver_structs.hpp"
 
 #include "dfpm/mesh.hpp"
 
 #include "reactmicp/systems/saturated_react/variables.hpp"
 #include "reactmicp/systems/saturated_react/init_variables.hpp"
 #include "reactmicp/systems/saturated_react/equilibrium_stagger.hpp"
 
 #include "dfpmsolver/parabolic_driver.hpp"
 #include "reactmicp/systems/saturated_react/transport_stagger.hpp"
 
 #include <iostream>
 #include "utils/log.hpp"
 
 TEST_CASE("Equilibrium stagger", "[SaturatedReact, chemistry, equilibrium, stagger]") {
 
     specmicp::logger::ErrFile::stream() = &std::cerr;
     specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
 
     specmicp::units::UnitsSet the_units;
     the_units.mass = specmicp::units::MassUnit::kilogram;
     the_units.length = specmicp::units::LengthUnit::centimeter;
 
     specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::uniform_mesh1d(5, 0.1, 5.0);
 
 
     specmicp::RawDatabasePtr the_database = leaching_database();
 
     std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition;
 
     list_initial_composition.push_back(initial_leaching_pb(the_database, 0.005, the_units));
     list_initial_composition.push_back(initial_blank_leaching_pb(the_database, 0.005, the_units));
 
     std::vector<specmicp::index_t> list_fixed_nodes = {0, };
 
     std::vector<int> index_initial_state = {1, 0, 0, 0, 0};
 
     specmicp::reactmicp::systems::satdiff::SaturatedVariablesPtr variables =
             specmicp::reactmicp::systems::satdiff::init_variables(the_mesh, the_database, the_units,
                                                                   list_fixed_nodes,
                                                                   list_initial_composition,
                                                                   index_initial_state);
 
     for (specmicp::index_t node: the_mesh->range_nodes())
     {
         variables->porosity(node) = 0.2;
         variables->diffusion_coefficient(node) = 1.0e-8;
     }
 
     SECTION("Initialisation") {
 
         specmicp::AdimensionalSystemConstraints constraints;
         constraints.set_saturated_system();
-        constraints.charge_keeper = 1;
+        constraints.charge_keeper = the_database->get_id_component("HO[-]");
 
         specmicp::AdimensionalSystemSolverOptions options;
         options.solver_options.maxstep = 10.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.use_scaling = true;
         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.fvectol = 1e-10;
         options.solver_options.steptol = 1e-14;
-        options.system_options.non_ideality_tolerance = 1e-10;
+        options.system_options.non_ideality_tolerance = 1e-12;
         options.units_set = the_units;
 
         variables->predictor() = variables->displacement();
 
         specmicp::reactmicp::systems::satdiff::EquilibriumStagger stagger(the_mesh->nb_nodes(), constraints, options);
 
         stagger.initialize_timestep(10, variables);
         stagger.restart_timestep(variables);
-
-        REQUIRE(variables->displacement() == variables->predictor());
+        REQUIRE((variables->displacement() - variables->predictor()).norm() == Approx(0).epsilon(1e-4));
     }
 
     SECTION("Transport + equilibrium") {
 
         specmicp::AdimensionalSystemConstraints constraints;
         constraints.set_saturated_system();
         constraints.charge_keeper = 1;
 
         specmicp::AdimensionalSystemSolverOptions options;
         options.solver_options.maxstep = 10.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-9;
-        options.solver_options.fvectol = 1e-6;
+        options.solver_options.fvectol = 1e-10;
         options.solver_options.steptol = 1e-14;
         options.system_options.non_ideality_tolerance = 1e-10;
         options.units_set = the_units;
 
         variables->predictor() = variables->displacement();
 
         specmicp::reactmicp::systems::satdiff::EquilibriumStagger equilibrium_stagger(the_mesh->nb_nodes() ,constraints, options);
 
         specmicp::reactmicp::systems::satdiff::SaturatedTransportStagger transport_stagger(variables, list_fixed_nodes);
 
         transport_stagger.initialize_timestep(10, variables);
         equilibrium_stagger.initialize_timestep(10, variables);
 
         int retcode;
         retcode = (int) transport_stagger.restart_timestep(variables);
         REQUIRE(retcode > 0);
         retcode = (int) equilibrium_stagger.restart_timestep(variables);
         REQUIRE(retcode > 0);
 
         REQUIRE(variables->displacement() != variables->predictor());
 
     }
 }
diff --git a/tests/reactmicp/systems/saturated_react/speciation_system.cpp b/tests/reactmicp/systems/saturated_react/speciation_system.cpp
index 870b1f7..bb4c4ba 100644
--- a/tests/reactmicp/systems/saturated_react/speciation_system.cpp
+++ b/tests/reactmicp/systems/saturated_react/speciation_system.cpp
@@ -1,131 +1,131 @@
 #include "database/database.hpp"
 #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 <iostream>
 
 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[-]"},
                                                 });
     thedatabase.keep_only_components({"H2O", "H[+]", "Ca[2+]", "Si(OH)4"});
     thedatabase.remove_gas_phases();
     thedatabase.swap_components(swapping);
     return thedatabase.get_database();
 }
 
 
 specmicp::AdimensionalSystemSolution initial_leaching_pb(
         specmicp::RawDatabasePtr raw_data,
         specmicp::scalar_t mult,
         specmicp::units::UnitsSet the_units)
 {
     specmicp::Formulation formulation;
 
     specmicp::scalar_t m_c3s = mult*0.7;
     specmicp::scalar_t m_c2s = mult*0.3;
     specmicp::scalar_t wc = 0.5;
     specmicp::scalar_t m_water = wc*1.0e-3*(
                   m_c3s*(3*56.08+60.08)
                 + m_c2s*(2*56.06+60.08)
                 );
 
     formulation.mass_solution = m_water;
     formulation.amount_minerals = {
         {"C3S", m_c3s},
         {"C2S", m_c2s},
     };
 
 
     specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
     specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
-    constraints.charge_keeper = 1;
+    constraints.charge_keeper = raw_data->get_id_component("HO[-]");;
     constraints.set_saturated_system();
 
     specmicp::AdimensionalSystemSolverOptions options;
     options.solver_options.maxstep = 100.0;
     options.solver_options.max_iter = 100;
     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-6;
+    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-10;
+    options.system_options.non_ideality_tolerance = 1e-12;
     options.units_set = the_units;
 
     specmicp::Vector x;
 
     specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
     solver.initialise_variables(x, 0.5, -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);
 }
 
 
 specmicp::AdimensionalSystemSolution initial_blank_leaching_pb(
         specmicp::RawDatabasePtr raw_data,
         specmicp::scalar_t mult,
         specmicp::units::UnitsSet the_units
         )
 {
     specmicp::Formulation formulation;
 
     specmicp::scalar_t m_c3s = mult*0.7;
     specmicp::scalar_t m_c2s = mult*0.3;
     specmicp::scalar_t wc = 0.5;
     specmicp::scalar_t m_water = wc*1e-3*(
                   m_c3s*(3*56.08+60.08)
                 + m_c2s*(2*56.06+60.08)
                 );
 
     formulation.mass_solution = m_water;
     formulation.concentration_aqueous =  {
         {"HO[-]",  1e-9},
         {"Ca[2+]", 1e-9},
         {"SiO(OH)3[-]", 1e-9}
     };
 
 
     specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
     specmicp::AdimensionalSystemConstraints constraints(total_concentrations);
-    constraints.charge_keeper = 1;
+    constraints.charge_keeper = raw_data->get_id_component("HO[-]");
     constraints.set_saturated_system();
 
     specmicp::AdimensionalSystemSolverOptions options;
     options.solver_options.maxstep = 10.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-9;
-    options.solver_options.fvectol = 1e-6;
+    options.solver_options.fvectol = 1e-10;
     options.solver_options.steptol = 1e-14;
-    options.system_options.non_ideality_tolerance = 1e-10;
+    options.system_options.non_ideality_tolerance = 1e-12;
     options.units_set = the_units;
 
     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);
 }
diff --git a/tests/reactmicp/systems/saturated_react/transport_program.cpp b/tests/reactmicp/systems/saturated_react/transport_program.cpp
index f9f0474..c148db8 100644
--- a/tests/reactmicp/systems/saturated_react/transport_program.cpp
+++ b/tests/reactmicp/systems/saturated_react/transport_program.cpp
@@ -1,87 +1,87 @@
 #include "catch.hpp"
 
 #include "speciation_system.hpp"
 #include "specmicp/adimensional/adimensional_system_solution.hpp"
 #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
 
 #include "dfpm/mesh.hpp"
 
 #include "reactmicp/systems/saturated_react/variables.hpp"
 #include "reactmicp/systems/saturated_react/init_variables.hpp"
 #include "reactmicp/systems/saturated_react/transport_program.hpp"
 
 #include "dfpmsolver/parabolic_driver.hpp"
 #include "reactmicp/systems/saturated_react/transport_stagger.hpp"
 
 TEST_CASE("Saturated transport program", "[SaturatedReact, transport, program, stagger]") {
 
     specmicp::units::UnitsSet the_units;
     the_units.mass = specmicp::units::MassUnit::kilogram;
     the_units.length = specmicp::units::LengthUnit::centimeter;
 
     specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::uniform_mesh1d(5, 0.1, 5.0);
 
 
     specmicp::RawDatabasePtr the_database = leaching_database();
 
     std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition;
 
     list_initial_composition.push_back(initial_leaching_pb(the_database, 0.005, the_units));
     list_initial_composition.push_back(initial_blank_leaching_pb(the_database, 0.005, the_units));
 
     std::vector<specmicp::index_t> list_fixed_nodes = {0, };
 
     std::vector<int> index_initial_state = {1, 0, 0, 0, 0};
 
     specmicp::reactmicp::systems::satdiff::SaturatedVariablesPtr variables =
             specmicp::reactmicp::systems::satdiff::init_variables(the_mesh, the_database, the_units,
                                                                   list_fixed_nodes,
                                                                   list_initial_composition,
                                                                   index_initial_state);
 
     for (specmicp::index_t node: the_mesh->range_nodes())
     {
         variables->porosity(node) = 0.2;
         variables->diffusion_coefficient(node) = 1.0e-8;
     }
 
     SECTION("Initialisation") {
 
         specmicp::reactmicp::systems::satdiff::SaturatedDiffusion program(variables, list_fixed_nodes);
 
         specmicp::Vector residuals;
 
-        REQUIRE(program.get_ndf() == 2*4);
-        REQUIRE(program.get_tot_ndf() == 2*4*5);
-        REQUIRE(program.get_neq() == 4*4);
+        REQUIRE(program.get_ndf() == 2*5); // 2*nb_comp
+        REQUIRE(program.get_tot_ndf() == 2*5*5); // 2*nb_comp*nb_nodes
+        REQUIRE(program.get_neq() == 4*3); // nb_nodes*active_comp
 
         program.compute_residuals(variables->displacement(), variables->velocity(), residuals);
 
 
 
     }
 
     SECTION("Solving the problem manually") {
         specmicp::reactmicp::systems::satdiff::SaturatedDiffusion program(variables, list_fixed_nodes);
 
         specmicp::dfpmsolver::ParabolicDriver<specmicp::reactmicp::systems::satdiff::SaturatedDiffusion> solver(program);
 
         variables->predictor() = variables->displacement();
 
          int retcode =  (int) solver.solve_timestep(100, variables->displacement());
 
          REQUIRE(retcode > 0);
     }
 
 
     SECTION("Transport stagger") {
         specmicp::reactmicp::systems::satdiff::SaturatedTransportStagger stagger(variables, list_fixed_nodes);
 
         stagger.initialize_timestep(100, variables);
 
         specmicp::reactmicp::solver::StaggerReturnCode retcode =  stagger.restart_timestep(variables);
 
         REQUIRE((int) retcode > 0);
 
     }
 }