diff --git a/examples/reactmicp/saturated_react/momas_benchmark.cpp b/examples/reactmicp/saturated_react/momas_benchmark.cpp
index 2e7e89d..34dd7d4 100644
--- a/examples/reactmicp/saturated_react/momas_benchmark.cpp
+++ b/examples/reactmicp/saturated_react/momas_benchmark.cpp
@@ -1,650 +1,677 @@
 #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 <iostream>
 #include <fstream>
 
 #include <chrono>
 #include <ctime>
 
 #include "database/selector.hpp"
 
 #include "reactmicp/solver/timestepper.hpp"
 
 #include "utils/io/reactive_transport.hpp"
 #include "utils/timer.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 = 0.25*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;
 
     AdimensionalSystemSolution solution = solver.get_raw_solution(variables);
     AdimensionalSystemSolutionExtractor sol(solution,
                                             raw_data,
                                             options.units_set);
 
     std::cout << sol.total_aqueous_concentration(2) << " - " << sol.total_immobile_concentration(2) << std::endl;
     std::cout << sol.density_water() << " - " << sol.volume_fraction_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; node<true_var->nb_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->diffusion_coefficient(node) = 330e-3;
+                //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->diffusion_coefficient(node) = 55e-3;
+                //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()
 {
     Timer global_timer;
     global_timer.start();
 
     std::ofstream output_gen("out_momas_out.txt");
 
     io::print_reactmicp_header(&output_gen);
 
     specmicp::logger::ErrFile::stream() = &output_gen; //&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 = 0.25*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 = 0.5*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-12;
     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 dx = 0.005;
     scalar_t tot_length = 2.1+dx;
 
     index_t nb_nodes = tot_length/dx +2;
 
     mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d(nb_nodes, dx, 1.0);
 
     std::vector<specmicp::AdimensionalSystemConstraints> list_constraints = {constraints_a,
                                                                              constraints_b,
                                                                              constraints_bc};
     std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition = {mat_a, mat_b, mat_bc};
     std::vector<int> index_initial_state(nb_nodes, 0);
     std::vector<int> 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<end_mat_b; ++node)
     {
         std::cout << "node b : " << node  << std::endl;
         index_initial_state[node] = 1;
         index_constraints[node] = 1;
     }
     index_initial_state[0] = 2;
     index_constraints[0] = 2;
 
     std::vector<specmicp::index_t> list_fixed_nodes = {0, };
     std::map<index_t, index_t> list_slaves_nodes = {{nb_nodes-1, nb_nodes-2}};
     std::vector<index_t> 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<reactmicp::systems::satdiff::EquilibriumStagger>(list_constraints, index_constraints, options);
 
     std::shared_ptr<reactmicp::systems::satdiff::SaturatedTransportStagger> transport_stagger =
             std::make_shared<reactmicp::systems::satdiff::SaturatedTransportStagger>(
                 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-14;
+    transport_options.step_tolerance = 1e-20;
+    transport_options.absolute_tolerance = 1e-6;
+
+
     transport_options.alpha = 1.0;
     //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<MoMasUspcaling>(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().good_enough_tolerance = 1.0;
 
-    solver.get_options().absolute_residuals_tolerance = 1e-12;
+    solver.get_options().absolute_residuals_tolerance = 1e-6;
     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;
 
 
 
     // mesh output
 
     std::ofstream output_mesh;
     output_mesh.open("out_momas_mesh.dat");
 
     output_mesh << "Node\tPosition" << std::endl;
     for (index_t node: the_mesh->range_nodes())
     {
         output_mesh << node << "\t" << the_mesh->get_position(node)-dx/2.0 << std::endl;
     }
     output_mesh.close();
 
 
     //std::cout << variables->displacement() << std::endl;
     std::ofstream output_iter("out_momas_iter.dat");
     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;
 
-    scalar_t dt=0.001;
-    reactmicp::solver::Timestepper timestepper(1e-3, 10.0, 100, 2);
+    scalar_t dt=0.0001;
+    reactmicp::solver::Timestepper timestepper(1e-4, 10.0, 5000, 2);
+    timestepper.get_options().alpha_average = 0.3;
+    timestepper.get_options().decrease_failure = 0.1;
+    timestepper.get_options().increase_factor = 1.05;
+    timestepper.get_options().decrease_factor = 0.50;
+    timestepper.get_options().iteration_lower_target = 1.001;
 
     io::print_reactmicp_performance_long_header(&output_iter);
 
     while (timestepper.get_total() < timestepper.get_total_target())
     {
         output_gen << "dt : " << dt << std::endl;
 
         Timer step_timer;
         step_timer.start();
         reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables);
         step_timer.stop();
 
-        if (retcode > reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
-        {
+       // if (retcode > reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
+        //{
             io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total()+dt, solver.get_perfs());
-        }
+        //}
         dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
         if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
         {
-            dt = 0.01;
+            dt = 0.0001;
             variables->reset_main_variables();
             retcode = solver.solve_timestep(dt, variables);
             dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
+            io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total(), solver.get_perfs());
         }
 
         output_gen  << "Time : " << timestepper.get_total() << std::endl;
         io::print_reactmicp_performance_short(&output_gen, solver.get_perfs());
 
         scalar_t total = timestepper.get_total();
 
-        if (cnt % 2== 0)
+        if (cnt % 10== 0)
         {
             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;
     }
 
-    reactmicp::solver::ReactiveTransportTimer timer = solver.get_timer();
-
-    io::print_reactmicp_timer(&output_gen, timer);
-
-    global_timer.stop();
-    output_gen << "Total computation time : " << global_timer.elapsed_time() << "s." << std::endl;
-
-    return EXIT_SUCCESS;
-
     variables->displacement()(1) = 0;
     variables->displacement()(2) = -2.0;
     variables->displacement()(3) = 0;
     variables->displacement()(4) = 2.0;
     variables->displacement()(6) = 0.0;
     variables->displacement()(7) = 0.0;
     variables->displacement()(8) = 0.0;
 
     timestepper.set_total_target(6000);
-    dt = 0.1;
+    dt = 1e-4;
+
+    transport_options.absolute_tolerance = 1e-10;
+    solver.get_options().absolute_residuals_tolerance = 1e-6;
 
     while (timestepper.get_total() < timestepper.get_total_target())
     {
-        output_gen << "dt : " << dt << std::endl;
-        reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables);
 
+        output_gen << "dt : " << dt << std::endl;
 
+        Timer step_timer;
+        step_timer.start();
+        reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables);
+        step_timer.stop();
 
+       // if (retcode > reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
+        //{
+            io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total()+dt, solver.get_perfs());
+        //}
         dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
+        if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
+        {
+            dt = 0.0001;
+            variables->reset_main_variables();
+            retcode = solver.solve_timestep(dt, variables);
+            dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
+            io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total(), solver.get_perfs());
+        }
 
-        output_gen  << "Time : " << timestepper.get_total()
-                   << "\nIter : " << solver.get_perfs().nb_iterations
-                   << "\nResiduals : " << solver.get_perfs().residuals << std::endl;
+        output_gen  << "Time : " << timestepper.get_total() << std::endl;
+        io::print_reactmicp_performance_short(&output_gen, solver.get_perfs());
 
         scalar_t total = timestepper.get_total();
 
-        if (cnt % 2== 0)
+        if (cnt % 10== 0)
         {
             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;
     }
 
+    reactmicp::solver::ReactiveTransportTimer timer = solver.get_timer();
+
+    io::print_reactmicp_timer(&output_gen, timer);
+
+    global_timer.stop();
+    output_gen << "Total computation time : " << global_timer.elapsed_time() << "s." << std::endl;
+
+
     return EXIT_SUCCESS;
 }