diff --git a/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp b/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp
index fa0f533..cb19e8f 100644
--- a/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp
+++ b/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp
@@ -1,118 +1,118 @@
 #include "equilibrium_stagger.hpp"
 #include "variables.hpp"
 
 #include "specmicp/adimensional/adimensional_system_solver.hpp"
 #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
 #include "utils/log.hpp"
 #include "reactmicp/solver/staggers_base/stagger_structs.hpp"
 
 #ifdef SPECMICP_USE_OPENMP
 #include <omp.h>
 #endif // SPECMICP_USE_OPENMP
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 // EquilibriumStagger::
 
 //! \brief Initialize the stagger at the beginning of an iteration
 void EquilibriumStagger::initialize_timestep(scalar_t dt, VariablesBasePtr var)
 {
     m_dt = dt;
 
     SaturatedVariablesPtr true_var = cast_var_from_base(var);
 
     // Initialize velocity using values from previous timestep
     for (index_t node=0; node<true_var->nb_nodes(); ++node)
     {
         if (true_var->is_fixed_composition(node)) continue;
         scalar_t alpha = 1.0;
         for (index_t component; component<true_var->nb_component(); ++component)
         {
             alpha = std::max(alpha,
                              dt*true_var->solid_concentration(node, component, true_var->chemistry_rate())
                                /true_var->solid_concentration(node, component, true_var->displacement())
                              );
         }
         auto solid_velocity = true_var->velocity().segment(
                     true_var->offset_node(node)+true_var->offset_solid_concentration(),
                     true_var->nb_component());
         auto solid_chemistry_rate = true_var->chemistry_rate().segment(
                     true_var->offset_node(node)+true_var->offset_solid_concentration(),
                     true_var->nb_component());
         solid_velocity = 1/alpha * solid_chemistry_rate;
     }
 }
 
 //! \brief Solve the equation for the timestep
 solver::StaggerReturnCode EquilibriumStagger::restart_timestep(VariablesBasePtr var)
 {
     SaturatedVariablesPtr true_var = cast_var_from_base(var);
 #ifdef SPECMICP_USE_OPENMP
 #pragma omp parallel
     {
 #pragma omp for schedule(runtime)
 #endif // SPECMICP_USE_OPENMP
-        for (index_t node=0; node<true_var->nb_component(); ++node)
+        for (index_t node=0; node<true_var->nb_nodes(); ++node)
         {
-            if (true_var->is_fixed_composition(false)) continue;
+            if (true_var->is_fixed_composition(node)) continue;
             solve_one_node(node, true_var);
         }
 #ifdef SPECMICP_USE_OPENMP
     }
 #endif // SPECMICP_USE_OPENMP
-    return solver::StaggerReturnCode::NotConvergedYet;
+    return solver::StaggerReturnCode::ResidualMinimized;
 }
 
 //!
 void EquilibriumStagger::solve_one_node(index_t node, SaturatedVariablesPtr true_var)
 {
     AdimensionalSystemConstraints constraints(m_constraints);
     constraints.total_concentrations = true_var->total_concentrations(node);
 
     AdimensionalSystemSolver adim_solver(true_var->get_database(),
                                          constraints,
                                          true_var->equilibrium_solution(node),
                                          m_options);
     Vector variables(true_var->equilibrium_solution(node).main_variables);
     micpsolver::MiCPPerformance perf =  adim_solver.solve(variables);
 
     micpsolver::MiCPSolverReturnCode retcode = perf.return_code;
     if (retcode <= micpsolver::MiCPSolverReturnCode::NotConvergedYet)
     {
         ERROR << "Failed to solve chemistry problem at node " << node
               << ", return code = " << static_cast<int>(retcode);
     }
 
     true_var->equilibrium_solution(node) = adim_solver.get_raw_solution(variables);
     AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node),
                                                   true_var->get_database(),
                                                   m_options.units_set);
 
     for (index_t component=0; component<true_var->nb_component(); ++component)
     {
         const scalar_t c_aq = extractor.density_water()*extractor.total_aqueous_concentration(component);
         true_var->aqueous_concentration(node, component, true_var->displacement()) = c_aq;
 
         const scalar_t c_aq_0 = true_var->aqueous_concentration(node, component, true_var->predictor());
         const scalar_t vel_aq = (c_aq - c_aq_0)/m_dt;
         true_var->aqueous_concentration(node, component, true_var->velocity()) = vel_aq;
 
         const scalar_t c_sol = extractor.total_solid_concentration(component);
         true_var->solid_concentration(node, component, true_var->displacement()) = c_sol;
 
         const scalar_t c_sol_0 = true_var->solid_concentration(node, component, true_var->predictor());
         const scalar_t vel_sol = (c_sol - c_sol_0)/m_dt;
         true_var->solid_concentration(node, component, true_var->velocity()) = vel_sol;
         true_var->solid_concentration(node, component, true_var->chemistry_rate()) =  vel_sol;
     }
 
 
 }
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
diff --git a/src/reactmicp/systems/saturated_react/react_solver.cpp b/src/reactmicp/systems/saturated_react/react_solver.cpp
index e69de29..8668cb2 100644
--- a/src/reactmicp/systems/saturated_react/react_solver.cpp
+++ b/src/reactmicp/systems/saturated_react/react_solver.cpp
@@ -0,0 +1,11 @@
+#include "react_solver.hpp"
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace satdiff {
+
+} // end namespace satdiff
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
diff --git a/src/reactmicp/systems/saturated_react/react_solver.hpp b/src/reactmicp/systems/saturated_react/react_solver.hpp
index e69de29..342c1f6 100644
--- a/src/reactmicp/systems/saturated_react/react_solver.hpp
+++ b/src/reactmicp/systems/saturated_react/react_solver.hpp
@@ -0,0 +1,18 @@
+#ifndef SPECMICP_REACTMICP_SYSTEMS_SATURATED_REACTSOLVER_HPP
+#define SPECMICP_REACTMICP_SYSTEMS_SATURATED_REACTSOLVER_HPP
+
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace satdiff {
+
+
+
+} // end namespace satdiff
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
+
+
+#endif // SPECMICP_REACTMICP_SYSTEMS_SATURATED_REACTSOLVER_HPP
diff --git a/src/reactmicp/systems/saturated_react/transport_program.cpp b/src/reactmicp/systems/saturated_react/transport_program.cpp
index 54a88d5..5e4ce0d 100644
--- a/src/reactmicp/systems/saturated_react/transport_program.cpp
+++ b/src/reactmicp/systems/saturated_react/transport_program.cpp
@@ -1,234 +1,241 @@
 #include "transport_program.hpp"
 #include "dfpm/meshes/mesh1d.hpp"
 #include "variables.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 //class SaturatedDiffusion::
 
 
 SaturatedDiffusion::SaturatedDiffusion(SaturatedVariablesPtr variables,
                    std::vector<index_t> list_fixed_nodes):
     m_mesh(variables->get_mesh()),
     m_variables(variables),
     m_is_in_residual_computation(false)
 {
     number_equations(list_fixed_nodes);
 }
 
 //! \brief Return the number of degrees of freedom per node
 index_t SaturatedDiffusion::get_ndf() const {return 2*m_variables->nb_component();}
 
 //! \brief Return the total number of degrees of freedom
 index_t SaturatedDiffusion::get_tot_ndf() const {return 2*m_variables->nb_component()*m_variables->get_mesh()->nb_nodes();}
 
 
 void SaturatedDiffusion::number_equations(std::vector<index_t> list_fixed_nodes)
 {
     m_ideq.resizeLike(m_variables->displacement());
     m_ideq.setZero();
     for (auto it=list_fixed_nodes.begin(); it<list_fixed_nodes.end(); ++it)
     {
-        for (index_t component=0; component<m_variables->nb_component(); ++component)
+        for (index_t component=1; component<m_variables->nb_component(); ++component)
         {
             m_ideq(m_variables->dof_aqueous_concentration(*it, component)) = no_equation;
         }
     }
 
     index_t neq = 0;
     for (index_t node=0; node<m_mesh->nb_nodes(); ++node)
     {
-        for (index_t component=0; component<m_variables->nb_component(); ++component)
+        m_ideq(m_variables->dof_aqueous_concentration(node, 0)) = no_equation;
+        for (index_t component=1; component<m_variables->nb_component(); ++component)
         {
             index_t dof = m_variables->dof_aqueous_concentration(node, component);
             if (m_ideq(dof) != no_equation)
             {
                 m_ideq(dof) = neq;
                 ++neq;
             }
             dof = m_variables->dof_solid_concentration(node, component);
             m_ideq(dof) = no_equation;
         }
     }
     m_neq = neq;
 
 }
 
 void SaturatedDiffusion::compute_residuals(
         const Vector& displacement,
         const Vector& velocity,
         Vector& residual)
 {
     residual = Vector::Zero(get_neq());
     m_is_in_residual_computation = true;
     for (index_t element: m_mesh->range_elements())
     {
-        for (index_t component=0; component<m_variables->nb_component(); ++component)
+        for (index_t component=1; component<m_variables->nb_component(); ++component)
         {
             Eigen::Vector2d element_residual;
             element_residual.setZero();
             residuals_element_component(element, component, displacement, velocity, element_residual);
             for (index_t en=0; en<2; ++en)
             {
                 const index_t node = m_mesh->get_node(element, en);
                 const index_t id = m_ideq(m_variables->dof_aqueous_concentration(node, component));
                 if (id != no_equation) {residual(id) += element_residual(en);}
             }
         }
     }
     m_is_in_residual_computation = false;
 }
 
 void SaturatedDiffusion::residuals_element_component(
         index_t element,
         index_t component,
         const Vector& displacement,
         const Vector& velocity,
         Eigen::Vector2d& element_residual
         )
 {
-
     const scalar_t mass_coeff_0 = -m_mesh->get_volume_cell_element(element, 0);
     const scalar_t mass_coeff_1 = -m_mesh->get_volume_cell_element(element, 1);
 
     const index_t node_0 = m_mesh->get_node(element, 0);
     const index_t node_1 = m_mesh->get_node(element, 1);
 
     const scalar_t diff_coeff = 1.0/(0.5/m_variables->diffusion_coefficient(node_0) +
                                      0.5/m_variables->diffusion_coefficient(node_1));
 
 
     scalar_t flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element)
                            * diff_coeff
                            ) ;
 
     const index_t dof_0 = m_variables->dof_aqueous_concentration(node_0, component);
     const index_t dof_1 = m_variables->dof_aqueous_concentration(node_1, component);
 
 
     scalar_t flux = flux_coeff*(displacement(dof_0) - displacement(dof_1));
 
     element_residual(0) = flux;
     element_residual(1) = - flux;
 
     if (m_is_in_residual_computation)
     {
         m_variables->aqueous_concentration(node_0, component,
                            m_variables->transport_rate()) += element_residual(0);
         m_variables->aqueous_concentration(node_1, component,
                            m_variables->transport_rate()) += element_residual(1);
     }
     // velocity
     element_residual(0) += mass_coeff_0*(velocity(dof_0)*m_variables->porosity(node_0)
                                          +m_variables->vel_porosity(node_0)*displacement(dof_0));
     element_residual(1) += mass_coeff_1*(velocity(dof_1)*m_variables->porosity(node_1)
                                          + m_variables->vel_porosity(node_1)*displacement(dof_1));
     // external rate
-    for (index_t en=0; en<2; ++en)
-    {
-        element_residual(en) += m_variables->aqueous_concentration(
-                    m_mesh->get_node(element, en), component,
-                    m_variables->chemistry_rate());
-    }
+    element_residual(0) +=  mass_coeff_0*m_variables->solid_concentration(node_0, component,
+                m_variables->chemistry_rate());
+    element_residual(1) +=  mass_coeff_1*m_variables->solid_concentration(node_1, component,
+                m_variables->chemistry_rate());
+
+
+//    for (index_t en=0; en<2; ++en)
+//    {
+//        element_residual(en) +=  mass_coeff_0*m_variables->solid_concentration(
+//                    m_mesh->get_node(element, en), component,
+//                    m_variables->chemistry_rate());
+//    }
 }
 
 
 void SaturatedDiffusion::compute_jacobian(
         Vector& displacement,
         Vector& velocity,
         Eigen::SparseMatrix<scalar_t>& jacobian,
         scalar_t alphadt
         )
 {
     dfpm::list_triplet_t jacob;
     const index_t ncomp = m_variables->nb_component();
     const index_t estimation = m_mesh->nb_nodes()*(ncomp*m_mesh->nen);
     jacob.reserve(estimation);
     for (index_t element: m_mesh->range_elements())
     {
         jacobian_element(element, displacement, velocity, jacob, alphadt);
     }
     jacobian = Eigen::SparseMatrix<scalar_t>(get_neq(), get_neq());
     jacobian.setFromTriplets(jacob.begin(), jacob.end());
 }
 
 void SaturatedDiffusion::jacobian_element(
         index_t element,
         Vector& displacement,
         Vector& velocity,
         dfpm::list_triplet_t& jacobian,
         scalar_t alphadt)
 {
-    for (index_t component=0; component<m_variables->nb_component(); ++component)
+    for (index_t component=1; component<m_variables->nb_component(); ++component)
     {
         Eigen::Vector2d element_residual_orig(Eigen::Vector2d::Zero());
         residuals_element_component(element, component, displacement, velocity, element_residual_orig);
 
         for (index_t en=0; en<2; ++en)
         {
             Eigen::Vector2d element_residual(Eigen::Vector2d::Zero());
 
             const index_t node = m_mesh->get_node(element, en);
             const index_t dof = m_variables->dof_aqueous_concentration(node, component);
             const index_t idc = m_ideq(dof);
 
             if (idc == no_equation) continue;
 
             const scalar_t tmp_v = velocity(dof);
             const scalar_t tmp_d = displacement(dof);
 
             scalar_t h = eps_jacobian*std::abs(tmp_v);
             if (h < 1e-4*eps_jacobian) h = eps_jacobian;
             velocity(dof) = tmp_v + h;
             h = velocity(dof) - tmp_v;
 
             displacement(dof) = tmp_d + alphadt*h;
 
             residuals_element_component(element, component, displacement, velocity, element_residual);
             velocity(dof) = tmp_v;
             displacement(dof) = tmp_d;
 
             for (index_t enr=0; enr<2; ++enr)
             {
                 const index_t noder = m_mesh->get_node(element, enr);
                 const index_t idr = m_ideq(m_variables->dof_aqueous_concentration(noder, component));
 
                 if (idr == no_equation) continue;
                 jacobian.push_back(dfpm::triplet_t(
                               idr,
                               idc,
                               (element_residual(enr) - element_residual_orig(enr))/h
                               ));
             }
         }
     }
 }
 
 //! \brief Update the solutions
 void SaturatedDiffusion::update_solution(const Vector& update,
                      scalar_t lambda,
                      scalar_t alpha_dt,
                      Vector& predictor,
                      Vector& displacement,
                      Vector& velocity)
 {
     for (index_t node: m_mesh->range_nodes())
     {
-        for (index_t component=0; component< m_variables->nb_component(); ++component)
+        for (index_t component=1; component< m_variables->nb_component(); ++component)
         {
             const index_t dof = m_variables->dof_aqueous_concentration(node, component);
             const index_t id = m_ideq(dof);
             if (id == no_equation) continue;
             velocity(dof) += lambda*update(id);
         }
     }
-    displacement = m_variables->predictor() + alpha_dt*velocity;
+    //displacement = m_variables->predictor() + alpha_dt*velocity;
+    displacement = predictor + alpha_dt*velocity;
 }
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
diff --git a/src/reactmicp/systems/saturated_react/transport_stagger.cpp b/src/reactmicp/systems/saturated_react/transport_stagger.cpp
index ebdeeb3..29c9af8 100644
--- a/src/reactmicp/systems/saturated_react/transport_stagger.cpp
+++ b/src/reactmicp/systems/saturated_react/transport_stagger.cpp
@@ -1,105 +1,114 @@
 #include "transport_stagger.hpp"
 #include "variables.hpp"
 #include "reactmicp/solver/staggers_base/stagger_structs.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 // class SaturatedTransportStagger::
 
 namespace internal {
 
 //! \brief Translate a ParabolicDriverReturnCode to a StaggerReturnCode
 solver::StaggerReturnCode Parabolic2StaggerReturnCode(dfpmsolver::ParabolicDriverReturnCode retcode)
 {
     using ParRC = dfpmsolver::ParabolicDriverReturnCode;
     using StaRC = solver::StaggerReturnCode;
 
     switch (retcode)
     {
     case ParRC::ResidualMinimized:
             return StaRC::ResidualMinimized;
             break;
     case ParRC::ErrorMinimized:
             return StaRC::ErrorMinimized;
             break;
     default: // error
         switch(retcode)
         {
         case ParRC::MaxIterations:
             return StaRC::MaximumIterationsReached;
             break;
         case ParRC::StationaryPoint:
             return StaRC::StationaryPoint;
             break;
         case ParRC::NotConvergedYet:
             return StaRC::NotConvergedYet;
             break;
         default:
             return StaRC::UnknownError;
         }
     }
 }
 
 } // end namespace internal
 
 SaturatedTransportStagger::SaturatedTransportStagger(SaturatedVariablesPtr variables,
                                                      std::vector<index_t> list_fixed_nodes):
     m_program(variables, list_fixed_nodes),
     m_solver(m_program)
 {
 
 }
 
 //! \brief Initialize the stagger at the beginning of an iteration
 void SaturatedTransportStagger::initialize_timestep(scalar_t dt, VariablesBasePtr var)
 {
     m_dt = dt;
 
     SaturatedVariablesPtr true_var = cast_var_from_base(var);
     true_var->predictor() = true_var->displacement();
 
     true_var->velocity().setZero();
     true_var->transport_rate().setZero();
     m_solver.initialize_timestep(dt, true_var->displacement());
+
+    Vector tmp = true_var->chemistry_rate();
+    true_var->chemistry_rate().setZero();
+    Eigen::VectorXd residuals;
+    m_program.compute_residuals(true_var->displacement(), true_var->velocity(), residuals);
+    m_residual_0 = residuals.norm();
+
+    true_var->chemistry_rate() = tmp;
+
     m_program.set_variables(true_var);
 }
 
 
 //! \brief Solve the equation for the timestep
 solver::StaggerReturnCode SaturatedTransportStagger::restart_timestep(VariablesBasePtr var)
 {
     SaturatedVariablesPtr true_var = cast_var_from_base(var);
     m_solver.velocity() = true_var->velocity();
     dfpmsolver::ParabolicDriverReturnCode retcode = m_solver.restart_timestep(true_var->displacement());
     // copy variables if successful
     if (retcode > dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet)
     {
         true_var->velocity() = m_solver.velocity();
     }
 
     return internal::Parabolic2StaggerReturnCode(retcode);
 }
 
 scalar_t SaturatedTransportStagger::get_update(VariablesBasePtr var)
 {
     return cast_var_from_base(var)->velocity().norm();
 }
 
 
 //! \brief Compute the residuals norm
 scalar_t SaturatedTransportStagger::get_residual(VariablesBasePtr var)
 {
     SaturatedVariablesPtr true_var = cast_var_from_base(var);
 
     Eigen::VectorXd residuals;
     m_program.compute_residuals(true_var->displacement(), true_var->velocity(), residuals);
     return residuals.norm();
 }
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
diff --git a/src/reactmicp/systems/saturated_react/transport_stagger.hpp b/src/reactmicp/systems/saturated_react/transport_stagger.hpp
index f62bf61..4b3e2c8 100644
--- a/src/reactmicp/systems/saturated_react/transport_stagger.hpp
+++ b/src/reactmicp/systems/saturated_react/transport_stagger.hpp
@@ -1,58 +1,57 @@
 #ifndef SPECMICP_REACTMICP_SYSTEMS_SATURATED_TRANSPORTSTAGGER_HPP
 #define SPECMICP_REACTMICP_SYSTEMS_SATURATED_TRANSPORTSTAGGER_HPP
 
 #include "variablesfwd.hpp"
 #include "reactmicp/solver/staggers_base/transport_stagger_base.hpp"
 
 #include "transport_program.hpp"
 #include "dfpmsolver/parabolic_driver.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 using VariablesBasePtr = solver::VariablesBasePtr;
 
 class SaturatedTransportStagger: public solver::TransportStaggerBase
 {
 public:
 
     SaturatedTransportStagger(SaturatedVariablesPtr variables,
                               std::vector<index_t> list_fixed_nodes);
 
     //! \brief Return the options of the solver
-    dfpmsolver::ParabolicDriverOptions options_solver() {return m_solver.get_options();}
+    dfpmsolver::ParabolicDriverOptions& options_solver() {return m_solver.get_options();}
 
     //! \brief Initialize the stagger at the beginning of an iteration
     void initialize_timestep(scalar_t dt, VariablesBasePtr var) override;
 
     //! \brief Solve the equation for the timestep
     solver::StaggerReturnCode restart_timestep(VariablesBasePtr var) override;
 
     //! \brief Compute the residuals norm
     scalar_t get_residual(VariablesBasePtr var) override;
     //! \brief Compute the residuals norm
     scalar_t get_residual_0(VariablesBasePtr var) override {
-        // ther is no difference at this stage,
         // TODO : a slightly more optimized version would use the value from the solver
-        return get_residual(var);
+        return m_residual_0;
     }
 
     //! \brief Obtain the update
     scalar_t get_update(VariablesBasePtr var) override;
 
 private:
     scalar_t m_dt;
-
+    scalar_t m_residual_0;
     SaturatedDiffusion m_program;
     dfpmsolver::ParabolicDriver<SaturatedDiffusion> m_solver;
 
 };
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
 
 #endif // SPECMICP_REACTMICP_SYSTEMS_SATURATED_TRANSPORTSTAGGER_HPP
diff --git a/src/reactmicp/systems/saturated_react/variables.cpp b/src/reactmicp/systems/saturated_react/variables.cpp
index 17d9cc1..8a1fcc6 100644
--- a/src/reactmicp/systems/saturated_react/variables.cpp
+++ b/src/reactmicp/systems/saturated_react/variables.cpp
@@ -1,27 +1,26 @@
 #include "variables.hpp"
 
 #include "dfpm/meshes/mesh1d.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 SaturatedVariables::SaturatedVariables(mesh::Mesh1DPtr the_mesh,
                    RawDatabasePtr the_database):
     m_mesh(the_mesh),
     m_database(the_database)
 {}
 
 Vector SaturatedVariables::total_concentrations(index_t node) {
-    index_t offset = offset_node(node);
 
-    return porosity(node)*displacement().segment(offset+offset_aqueous_concentration(), nb_component())
-            + displacement().segment(offset+offset_solid_concentration(), nb_component());
+    return porosity(node)*displacement().segment(offset_aqueous_concentration(node), nb_component())
+            + displacement().segment(offset_solid_concentration(node), nb_component());
 }
 
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
diff --git a/src/reactmicp/systems/saturated_react/variables.hpp b/src/reactmicp/systems/saturated_react/variables.hpp
index fc2a21b..cdf7a57 100644
--- a/src/reactmicp/systems/saturated_react/variables.hpp
+++ b/src/reactmicp/systems/saturated_react/variables.hpp
@@ -1,198 +1,201 @@
 #ifndef SPECMICP_REACTMICP_SYSTEMS_SATURATED_VARIABLES_HPP
 #define SPECMICP_REACTMICP_SYSTEMS_SATURATED_VARIABLES_HPP
 
 #include "database.hpp"
 #include "reactmicp/solver/staggers_base/variables_base.hpp"
 #include "specmicp/adimensional/adimensional_system_solution.hpp"
 #include <vector>
 
 // forward declaration
 // ===================
 
 #include "dfpm/meshes/mesh1dfwd.hpp"
 namespace specmicp {
 namespace reactmicp {
 namespace solver {
     using VariablesBasePtr = std::shared_ptr<VariablesBase>;
 }
 namespace systems {
 namespace satdiff {
 class SaturatedVariablesFactory;
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
 
 
 // Class declaration
 // =================
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 //! \brief Variables for the saturated reactive transport system
 //!
 //! Contain all the variables that need to be shared between the staggers
 class SaturatedVariables: public solver::VariablesBase
 {
     // SaturatedVariablesFactory should be the class to use to inialize
     // the variables correctly
     friend class SaturatedVariablesFactory;
 
 public:
     SaturatedVariables(mesh::Mesh1DPtr the_mesh,
                        RawDatabasePtr the_database);
 
     //! \brief Return the mesh
     mesh::Mesh1DPtr get_mesh() {return m_mesh;}
     //! \brief Return the database
     RawDatabasePtr get_database() {return m_database;}
 
     //! \brief Return the number of components
     index_t nb_component() {return m_database->nb_component;}
     //! \brief Return the number of nodes
     index_t nb_nodes() {return m_is_fixed_composition.size();}
     //! \brief Return true if 'node' has a fixed composition
     index_t is_fixed_composition(index_t node) {return m_is_fixed_composition[node];}
 
     // Main variables
     // ==============
 
     //! \brief Return the main variable vector
     Vector& displacement() {return m_displacement;}
     //! \brief Return the main variable vector at the beginning of the timestep
     Vector& predictor() {return m_predictor;}
     //! \brief Return the velocity of the main variables
     Vector& velocity() {return m_velocity;}
     //! \brief Return the rate of change of the main variables due to the transport operator
     Vector& transport_rate() {return m_transport_rate;}
     //! \brief Return the rate of change of the main variables due to the chemistry operator
     Vector& chemistry_rate() {return m_chemistry_rate;}
 
     // Access to main variables
     // ========================
 
     //! \brief Return the number of degree of freedom (per node) in the main variables vector
     index_t ndf() {return 2*m_database->nb_component;}
     //! \brief Return the offset of 'node' in the main variables vector
     index_t offset_node(index_t node) {return node*ndf();}
     //! \brief Return the offset of the aqueous concentration variables in the main variables vector
     index_t offset_aqueous_concentration() {return 0;}
     //! \brief Return the offset of the aqueous concentrations variables in the main variables vector
     index_t offset_aqueous_concentration(index_t node) {
         return offset_aqueous_concentration()+offset_node(node);}
     //! \brief Return the offset of the solid concentration variables in the main variables vector
     index_t offset_solid_concentration() {return m_database->nb_component;}
     //! \brief Return the offset of the solid concentrations variables in the main variables vector
     index_t offset_solid_concentration(index_t node) {
         return offset_solid_concentration()+offset_node(node);}
 
     //! \brief Return the degree of freedom number for the aqueous concentration of 'component' at 'node'
     index_t dof_aqueous_concentration(index_t node, index_t component) {
         return (component + offset_aqueous_concentration(node));
     }
     //! \brief Return the degree of freedom number for the solid concentration of 'component' at 'node'
     index_t dof_solid_concentration(index_t node, index_t component) {
         return (component + offset_solid_concentration(node));
     }
     //! \brief Return the aqueous concentration of 'component' at 'node' in 'var'
     //!
     //! 'var' is any of the main variables vector, it may be a velocity vector
     scalar_t& aqueous_concentration(index_t node, index_t component, Vector& var) {
         return var(dof_aqueous_concentration(node, component));
     }
     //! \brief Return the solid concentration of 'component' at 'node' in 'var'
     //!
     //! 'var' is any of the main variables vector, it may be a velocity vector
     scalar_t& solid_concentration(index_t node, index_t component, Vector& var){
         return var(dof_solid_concentration(node, component));
     }
     //! \brief Return a vector containing the total concentrations computed from the main variables
     //!
     //! This is to be used to restart the chemistry computation
     Vector total_concentrations(index_t node);
 
     // Equilibrium
     // ===========
 
     //! \brief Returh the solution of the speciation solver at 'node'
     AdimensionalSystemSolution& equilibrium_solution(index_t node) {
         return m_equilibrium_solutions[node];
     }
 
     // Upscaling
     // =========
     //! \brief Return the offset for 'node' in the upscaling variables vector
     index_t offset_node_upscaling(index_t node) {return ndf_upscaling()*node;}
 
     //! \brief Return the number fo degree of freedom (per node) for the upscaling vector
     index_t ndf_upscaling() {return 3;}
     //! \brief Return the degree of freedom for the porosity at 'node'
     index_t dof_porosity(index_t node) {return 0 + offset_node_upscaling(node);}
     //! \brief Return the degree of freedom of the porosity velocity at 'node'
     index_t dof_vel_porosity(index_t node) {return 1 + offset_node_upscaling(node);}
     //! \brief Return the degree of freedom of the diffusion coefficient at 'node'
-    index_t dof_diffusion_coefficient(index_t node) {return 1 + offset_node_upscaling(node);}
+    index_t dof_diffusion_coefficient(index_t node) {return 2 + offset_node_upscaling(node);}
 
     //! \brief Return the porosity at 'node'
     scalar_t& porosity(index_t node) {return m_upscaling(dof_porosity(node));}
     //! \brief Return the rate of change of the porosity at 'node'
     scalar_t& vel_porosity(index_t node) {return m_upscaling(dof_vel_porosity(node));}
     //! \brief Return the diffusion coefficient at 'node'
     scalar_t& diffusion_coefficient(index_t node) {return m_upscaling(dof_diffusion_coefficient(node));}
 
+    //! \brief Return the vector of upscaling variables
+    Vector& upscaling_variables() {return m_upscaling;}
+
 private:
 
     // ############ //
     //  Attributes  //
     // ############ //
 
     mesh::Mesh1DPtr m_mesh;
     RawDatabasePtr m_database;
 
     std::vector<bool> m_is_fixed_composition;
 
     // Main variables
     // ==============
 
     Vector m_displacement;
     Vector m_predictor;
     Vector m_velocity;
     Vector m_transport_rate;
     Vector m_chemistry_rate;
 
     // Equilibrium
     // ===========
     std::vector<AdimensionalSystemSolution> m_equilibrium_solutions;
 
     // Upscaling
     // =========
 
     Vector m_upscaling;
 };
 
 //! \brief typedef of a shared pointer of a SaturatedVariables
 using SaturatedVariablesPtr = std::shared_ptr<SaturatedVariables>;
 
 // Casting function
 // =================
 
 //! \brief Static cast to a SaturatedVariablesPtr
 inline SaturatedVariablesPtr cast_var_from_base(solver::VariablesBasePtr var)
 {
     return std::static_pointer_cast<SaturatedVariables>(var);
 }
 //! \brief Static cast from a SaturatedVariablesPtr
 inline solver::VariablesBasePtr cast_var_to_base(SaturatedVariablesPtr var)
 {
     return std::static_pointer_cast<solver::VariablesBase>(var);
 }
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
 
 #endif // SPECMICP_REACTMICP_SYSTEMS_SATURATED_VARIABLES_HPP