diff --git a/examples/reactmicp/saturated_react/react_leaching.cpp b/examples/reactmicp/saturated_react/react_leaching.cpp
index 07cc313..e5761eb 100644
--- a/examples/reactmicp/saturated_react/react_leaching.cpp
+++ b/examples/reactmicp/saturated_react/react_leaching.cpp
@@ -1,448 +1,451 @@
 
 #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 <iostream>
 #include <fstream>
 
 using namespace specmicp;
 
 // Initial conditions
 // ==================
 
 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[-]"},
                                               {"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::AdimensionalSystemSolver solver(raw_data, constraints, options);
     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}
     };
 
 
     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 = 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-8;
     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.5;
     x.segment(1, raw_data->nb_component-1).setConstant(-3.0);
 
     AdimensionalSystemSolver solver(raw_data, constraints, options);
     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.4*extractor.density_water()*extractor.total_aqueous_concentration(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;
         for (index_t node=0; node<true_var->nb_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; node<true_var->nb_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; node<true_var->nb_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) = 1e-7;//std::min(1e4*std::exp(9.85*porosity-29.08),2.2e-5);
+        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.0100;
+    scalar_t dx = 0.0025;
     index_t additional_nodes = 1;
 
     radius = radius + additional_nodes*dx;
 
-    index_t nb_nodes = 10;
+    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);
 
-    mesh::Mesh1DPtr the_mesh = mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, dx, height);
     std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition = {initial, blank};
     std::vector<int> index_initial_state(nb_nodes, 0);
     for (int i=0; i< additional_nodes; ++i)
         index_initial_state[i] = 1;
     std::vector<specmicp::index_t> 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 = 10.0;
+    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-8;
-    options.solver_options.steptol = 1e-12;
+    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<reactmicp::systems::satdiff::EquilibriumStagger>(constraints, options);
 
     std::shared_ptr<reactmicp::systems::satdiff::SaturatedTransportStagger> transport_stagger =
             std::make_shared<reactmicp::systems::satdiff::SaturatedTransportStagger>(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<LeachingUspcaling>(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 < 20*3600*24)
+    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 % 100 == 0)
+        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;
     }
 }