diff --git a/src/reactmicp/meshes/mesh1d.hpp b/src/reactmicp/meshes/mesh1d.hpp
index 8d796b7..863919a 100644
--- a/src/reactmicp/meshes/mesh1d.hpp
+++ b/src/reactmicp/meshes/mesh1d.hpp
@@ -1,84 +1,86 @@
 /*-------------------------------------------------------
 
  - Module : reactmicp/meshes
  - File : mesh1d.hpp
  - Author : Fabien Georget
 
 Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
     * Neither the name of the Princeton University nor the
       names of its contributors may be used to endorse or promote products
       derived from this software without specific prior written permission.
 
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 ---------------------------------------------------------*/
 
 #ifndef SPECMICP_REACTMICP_MESH_MESH1D_HPP
 #define SPECMICP_REACTMICP_MESH_MESH1D_HPP
 
 #include "reactmicp/common.hpp"
 #include "boost/range/irange.hpp"
 #include <iostream>
 
 namespace specmicp {
 namespace reactmicp {
 namespace mesh {
 
 using range_t = boost::iterator_range<boost::range_detail::integer_iterator<ind_t>>;
 
 //! a one dimensional mesh
 class Mesh1D
 {
 public:
     Mesh1D(ind_t elements, double dx, double cross_section):
         m_dx(dx), m_crosssection(cross_section), m_nel(elements)
     {}
 
     const int nen = 2;
 
     ind_t nel() {return m_nel;}
     ind_t nnodes() {return m_nel+1;}
-    ind_t get_node(ind_t element, ind_t index) {return element+index;}
+    ind_t get_node(ind_t element, ind_t index) {
+        assert(index < nen);
+        return element+index;}
     double get_dx(ind_t element) {return m_dx;}
 
     double crosssection() { return m_crosssection;}
     double cell_volume(ind_t node) {
         if (node ==0 or node == nnodes()-1) return m_dx*m_crosssection/2;
         else return m_dx*m_crosssection;
     }
 
     range_t range_elements() {return boost::irange(0, nel());}
     range_t range_nodes() {return boost::irange(0, nnodes());}
     range_t range_nen() {return boost::irange(0, nen);}
 
 private:
 
     double m_dx;
     double m_crosssection;
     ind_t m_nel;
 
 };
 
 } // end namespace mesh
 } // end namespace reactmicp
 } // end namespace specmicp
 
 #endif // SPECMICP_REACTMICP_MESH_MESH1D_HPP
diff --git a/src/reactmicp/micpsolvers/micpsolver.hpp b/src/reactmicp/micpsolvers/micpsolver.hpp
index d3fe521..4ca2953 100644
--- a/src/reactmicp/micpsolvers/micpsolver.hpp
+++ b/src/reactmicp/micpsolvers/micpsolver.hpp
@@ -1,98 +1,106 @@
 #ifndef SPECMICP_REACTMICP_MICPSOLVER_MICPSOLVER_HPP
 #define SPECMICP_REACTMICP_MICPSOLVER_MICPSOLVER_HPP
 
 #include <memory>
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 
 #include "reactmicp/common.hpp"
 #include "micpsolver_structs.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace micpsolver  {
 
 //! MiCP solver for reactive transport
 template <class program_t>
 class MiCPParabolicSolver
 {
 public:
 
     using program_ptr_t = std::shared_ptr<program_t>;
 
     MiCPParabolicSolver(program_ptr_t the_program):
         m_program(the_program)
     {}
 
     //! \brief Solve the program
     MiCPParabolicSolverReturnCode solve_one_timestep(double dt,
                                                     ParabolicVariables& x);
 
 
     int linesearch(Eigen::VectorXd& update, ParabolicVariables& x,
                    const Eigen::VectorXd &predictor, double dt);
-    MiCPParabolicSolverReturnCode compute_search_direction(Eigen::VectorXd& update);
+    MiCPParabolicSolverReturnCode compute_search_direction(Eigen::VectorXd& update,
+                                                           ParabolicVariables &x,
+                                                           double dt);
 
     //! \brief Compute and reformulate the residuals
     void setup_residuals(const ParabolicVariables& x);
     //! \brief Compute the residuals
     void compute_residuals(const ParabolicVariables& x,
                            Eigen::VectorXd& residual);
     //! \brief Reformulate the residual
     void reformulate_residuals(const ParabolicVariables& x);
     //! \brief Reformulate the residual *in place*
     void reformulate_residuals_inplace(const ParabolicVariables& x, Eigen::VectorXd& residual);
     //! \brief Setup the jacobian
-    void setup_jacobian(ParabolicVariables& x, const double dt);
+    void setup_jacobian(ParabolicVariables& x,
+                        const double dt,
+                        Eigen::SparseMatrix<double, Eigen::RowMajor> &jacobian);
     //! \brief Compute the jacobian
-    void compute_jacobian(ParabolicVariables& x, const double dt);
+    void compute_jacobian(ParabolicVariables& x,
+                          const double dt,
+                           Eigen::SparseMatrix<double, Eigen::RowMajor> &jacobian);
     //! \brief Reformulate the jacobian
-    void reformulate_jacobian(const ParabolicVariables& x);
+    void reformulate_jacobian(const ParabolicVariables& x,
+                              Eigen::SparseMatrix<double, Eigen::RowMajor> &jacobian, double dt);
 
     MiCPParabolicSolverReturnCode check_convergence(int nb_iterations,
                                                 const Eigen::VectorXd& update,
                                                 const ParabolicVariables& solution,
                                                 const Eigen::VectorXd& residuals,
                                                 bool may_have_converged);
 
 
     //! \brief Return the program
     program_ptr_t get_program() {return m_program;}
     //! \brief Return a reference to the options
     MiCPParabolicSolverOptions& get_options() {return m_options;}
     //! \brief Return a reference to the performance
     MiCPParabolicSolverPerformance& get_performance() {return m_performance;}
 
     //! \brief Scale the step if it is too long
     double is_step_too_long(Eigen::VectorXd& update);
 
     //! \brief Return the number of equations
     double get_neq() {return m_program->get_neq();}
 private:
     program_ptr_t m_program;
     MiCPParabolicSolverOptions m_options;
     MiCPParabolicSolverPerformance m_performance;
     Eigen::VectorXd m_residuals;
     Eigen::VectorXd m_phi_residuals;
     Eigen::VectorXd m_gradient;
     Eigen::VectorXd m_predictor;
-    Eigen::SparseMatrix<double, Eigen::RowMajor> m_jacobian;
+    Eigen::VectorXd m_cp_variables;
 
     std::vector<double> m_max_merits;   //!< Contains the m best value of the merit function
     bool m_gradient_step_taken;
 
+    double m_norm_0;
 };
 
 } // end namespace micpsolver
 } // end namespace reactmicp
 } // end namespace specmicp
 
 // ============== //
 //                //
 // implementation //
 //                //
 // ============== //
 #include "micpsolver.inl"
 
 
 #endif // SPECMICP_REACTMICP_MICPSOLVER_MICPSOLVER_HPP
diff --git a/src/reactmicp/micpsolvers/micpsolver.inl b/src/reactmicp/micpsolvers/micpsolver.inl
index cb00718..aad6955 100644
--- a/src/reactmicp/micpsolvers/micpsolver.inl
+++ b/src/reactmicp/micpsolvers/micpsolver.inl
@@ -1,452 +1,536 @@
 #include "micpsolver.hpp" // for syntaxic coloration only
 
 #include "micpsolver/ncp_function.hpp"
 #include "utils/log.hpp"
 #include <iostream>
+    #include <unsupported/Eigen/SparseExtra>
 
 namespace specmicp {
 namespace reactmicp {
 namespace micpsolver  {
 
 // ============================= //
 //                               //
 //        Main algorithm         //
 //        --------------         //
 //                               //
 // ============================= //
 
 template <class program_t>
 MiCPParabolicSolverReturnCode MiCPParabolicSolver<program_t>::solve_one_timestep(double dt,
                                        ParabolicVariables& x)
 {
     int cnt = 0;
     m_residuals.resize(get_neq());
-    setup_residuals(x); // << is it necessary ??
     MiCPParabolicSolverReturnCode retcode = MiCPParabolicSolverReturnCode::NotConvergedYet;
     Eigen::VectorXd update(get_neq());
     Eigen::VectorXd predictor;
     m_program->set_predictor(x, predictor, get_options().alpha, dt);
+    setup_residuals(x);
+    m_norm_0 = m_phi_residuals.norm();
     while (retcode == MiCPParabolicSolverReturnCode::NotConvergedYet)
     {
+        //std::cout << "iteration : " << cnt << std::endl;
         DEBUG << "Iteration : " << cnt;
         // (S.0) hook
         const bool may_have_converged = get_program()->hook_start_iteration(x, m_phi_residuals.norm());
         // (S.1) compute the residuals
         setup_residuals(x);
+        //std::cout << "Residuals " << std::endl << m_phi_residuals << std::endl;
         // (S.2) check the convergence
         retcode = check_convergence(cnt, update, x, m_phi_residuals, may_have_converged);
         get_performance().return_code = retcode;
+        //std::cout << "Retcode : " << (int) retcode << std::endl;
+        if (retcode != MiCPParabolicSolverReturnCode::NotConvergedYet) break;
         ++cnt;
         // (S.3) compute the jacobian
-        setup_jacobian(x, dt);
-        compute_search_direction(update);
+        compute_search_direction(update, x, dt);
         // (S.4) linesearch and update
         linesearch(update, x, predictor, dt);
     }
     get_performance().nb_iterations = cnt;
     return retcode;
 }
 
 
 
 template <class program_t>
 MiCPParabolicSolverReturnCode MiCPParabolicSolver<program_t>::check_convergence(
                                             int nb_iterations,
                                             const Eigen::VectorXd& update,
                                             const ParabolicVariables& solution,
                                             const Eigen::VectorXd& residuals,
                                             bool may_have_converged)
 {
     MiCPParabolicSolverReturnCode termcode = MiCPParabolicSolverReturnCode::NotConvergedYet;
     const double norm_residuals = residuals.lpNorm<Eigen::Infinity>();
-    WARNING << "Residuals : " << nb_iterations << " - " << norm_residuals;
-    if (norm_residuals < get_options().residuals_tolerance)
+    WARNING << "Residuals : " << nb_iterations << " - " << norm_residuals/m_norm_0;
+    if (norm_residuals/m_norm_0 < get_options().residuals_tolerance)
     {
         if (may_have_converged == true)
             termcode = MiCPParabolicSolverReturnCode::ResidualMinimized;
     }
     else if (nb_iterations > 0 and update.lpNorm<Eigen::Infinity>() < get_options().step_tolerance)
     {
         if (norm_residuals > get_options().threshold_stationary_point)
         {
             ERROR << "Stationary point detected !";
             termcode = MiCPParabolicSolverReturnCode::StationaryPoint;
         }
         if (may_have_converged == true)
         {
             WARNING << "Error is minimized - may indicate a stationnary point";
             termcode = MiCPParabolicSolverReturnCode::ErrorMinimized;
         }
     }
     else if (nb_iterations >  get_options().maximum_iterations)
     {
         ERROR << "Maximum number of iteration reached (" << get_options().maximum_iterations << ")";
         termcode = MiCPParabolicSolverReturnCode::MaxIterations;
     }
     else if (get_performance().maximum_step_taken)
     {
         ++get_performance().nb_consecutive_maximum_step_taken;
         ++get_performance().nb_maximum_step_taken;
         if (get_performance().nb_consecutive_maximum_step_taken == get_options().maximum_iterations_at_maximum_step) {
             ERROR << "Divergence detected - Maximum step length taken two many times";
             termcode = MiCPParabolicSolverReturnCode::MaxStepTakenTooManyTimes;
         }
     }
     else
     {
         get_performance().nb_consecutive_maximum_step_taken = 0;
     }
     return termcode;
 }
 
 // ============================== //
 //                                //
 //     Residuals and Jacobian     //
 //     ----------------------     //
 //                                //
 // ============================== //
 
 // Residuals
 // =========
 
 template <class program_t>
 void MiCPParabolicSolver<program_t>::setup_residuals(const ParabolicVariables& x)
 {
     compute_residuals(x, m_residuals);
     reformulate_residuals(x);
 }
 
 template <class program_t>
 void MiCPParabolicSolver<program_t>::compute_residuals(const ParabolicVariables& x,
                                                     Eigen::VectorXd& residuals)
 {
+    m_program->copy_cp_variables(x, m_cp_variables);
     get_program()->get_residuals(x, residuals);
     get_performance().nb_call_residuals += 1;
 }
 
 // reformulation
 // -------------
 
 template <class program_t>
 void MiCPParabolicSolver<program_t>::reformulate_residuals(const ParabolicVariables &x)
 {
     const std::vector<bool>& is_cp = m_program->micp_index();
     m_phi_residuals.resizeLike(m_residuals);
+    m_phi_residuals.setZero();
     for (int i=0; i<m_program->get_neq(); ++i)
     {
         if (is_cp[i])
         {
             m_phi_residuals(i) = specmicp::micpsolver::penalized_fisher_burmeister(
-                        x.displacement(i),
+                        m_cp_variables(i),
                         m_residuals(i),
                         get_options().penalization_factor
                         );
         }
+        else
+        {
+            m_phi_residuals(i) = m_residuals(i);
+        }
     }
 }
 
 // this is mainly used in the linesearch
 template <class program_t>
 void MiCPParabolicSolver<program_t>::reformulate_residuals_inplace(const ParabolicVariables& x,
                                                Eigen::VectorXd& residual)
 {
     const std::vector<bool>& is_cp = m_program->micp_index();
 
     for (int i=0; i<m_program->get_neq(); ++i)
     {
         if (is_cp[i])
         {
-            residual(i) = specmicp::micpsolver::penalized_fisher_burmeister(
-                        x.displacement(i),
+            residual.coeffRef(i) = specmicp::micpsolver::penalized_fisher_burmeister(
+                        m_cp_variables(i),
                         residual(i),
                         get_options().penalization_factor
                         );
         }
     }
 }
 
 // Jacobian
 // ========
 
 template <class program_t>
-void MiCPParabolicSolver<program_t>::setup_jacobian(ParabolicVariables& x, const double dt)
+void MiCPParabolicSolver<program_t>::setup_jacobian(
+        ParabolicVariables& x,
+        const double dt,
+        Eigen::SparseMatrix<double, Eigen::RowMajor>& jacobian)
 {
-    compute_jacobian(x, dt);
-    reformulate_jacobian(x);
-    m_jacobian.makeCompressed();
+    compute_jacobian(x, dt, jacobian);
+    reformulate_jacobian(x, jacobian, dt);
+    jacobian.makeCompressed();
 }
 
 template <class program_t>
-void MiCPParabolicSolver<program_t>::compute_jacobian(ParabolicVariables& x, const double dt)
+void MiCPParabolicSolver<program_t>::compute_jacobian(ParabolicVariables& x,
+                                                      const double dt,
+                                                      Eigen::SparseMatrix<double, Eigen::RowMajor> &jacobian)
 {
-    get_program()->get_jacobian(x, m_residuals, m_jacobian, get_options().alpha*dt);
+    get_program()->get_jacobian(x, m_residuals, jacobian, get_options().alpha*dt);
     get_performance().nb_call_jacobian += 1;
 }
 
 // Reformulation
 // -------------
 template <class program_t>
-void MiCPParabolicSolver<program_t>::reformulate_jacobian(const ParabolicVariables& x)
+void MiCPParabolicSolver<program_t>::reformulate_jacobian(const ParabolicVariables& x,
+                                                          Eigen::SparseMatrix<double, Eigen::RowMajor> &jacobian,
+                                                          double dt)
 {
     // see Facchinei and Pang for more details
 
     // adapted to sparse matrix in ROW major
     // i.e. the outer direction is a row (the algorithm proceed row by row)
 
     // Iterating over the non zero element :
     // see Eigen doc http://eigen.tuxfamily.org/dox/group__TutorialSparse.html
     // and https://forum.kde.org/viewtopic.php?f=74&t=122606&sid=a5b017a2e6dad7bab6f09e90680275b5
 
-    m_jacobian.makeCompressed(); // should be the case, but just to be sure
+    jacobian.makeCompressed(); // should be the case, but just to be sure
+    //std::cout << "residuals" << std::endl;
+    //std::cout << m_phi_residuals << std::endl;
+    //std::cout << "Jacobian " << std::endl;
+    //std::cout << jacobian.toDense() << std::endl;
     const std::vector<bool>&  is_cp = m_program->micp_index();
     const double lambda = get_options().penalization_factor;
     for (int i=0; i<get_neq(); ++i)
     {
         if (not is_cp[i]) continue;
-
         double c, d;
-        if (x.displacement(i) == 0 or m_residuals(i) == 0)
+        const double z = 1.0;
+        if (m_cp_variables(i) == 0 and m_residuals(i) == 0)
         {
             double gpdotz = 0; // grad(Fi)^t z
-            for (int k=0; k<m_jacobian.outerSize(); ++k)
-                for (SparseMatrix::InnerIterator it(m_jacobian,k); it; ++it)
+            for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(jacobian,i); it; ++it)
+            {
+                if (m_cp_variables(it.col()) == 0 and m_residuals(it.col()) == 0)
                 {
-                    if (x.displacement(it.col()) == 0 and m_residuals(it.col()) == 0)
-                    {
-                        gpdotz += it.value();
-                    }
+                    gpdotz += z*it.value();
                 }
-            double s = 1.0 + std::abs(gpdotz); // z(i) = 1.0
-            s = s * std::sqrt(1.0/(s*s) + gpdotz*gpdotz/(s*s));
-            c = lambda*(1.0/s - 1);
+            }
+            double s = z + std::abs(gpdotz);
+            s = s * std::sqrt(z*z/(s*s) + gpdotz*gpdotz/(s*s));
+            c = lambda*(z/s - 1);
             d = lambda*(gpdotz/s -1);
 
         }
         else
         {
-            double s = std::abs(x.displacement(i)) + std::abs( m_residuals(i));
-            s = s * std::sqrt(x.displacement(i)*x.displacement(i)/(s*s) +  m_residuals(i)* m_residuals(i)/(s*s));
-            c = lambda*(x.displacement(i)/s - 1);
-            d = lambda*( m_residuals(i)/s - 1);
-            if ((lambda <1) and ( m_residuals(i) > 0)  and (x.displacement(i) >0))
+            double s = std::abs(m_cp_variables(i)) + std::abs( m_residuals(i));
+            s = s * std::sqrt((m_cp_variables(i)*m_cp_variables(i))/(s*s) +  (m_residuals(i)*m_residuals(i))/(s*s));
+            c = lambda*(m_cp_variables(i)/s - 1.0);
+            d = lambda*( m_residuals(i)/s - 1.0);
+            if ((lambda < 1.0) and ( m_residuals(i) > 0.0)  and (m_cp_variables(i) > 0.0))
             {
-                c -= (1-lambda)* m_residuals(i);
-                d -= (1-lambda)*x.displacement(i);
+                c -= (1.0-lambda)* m_residuals(i);
+                d -= (1.0-lambda)*m_cp_variables(i);
             }
         }
+        c *= get_options().alpha*dt;
         // update the jacobian
-        for (int k=0; k<m_jacobian.outerSize(); ++k)
+        //std::cout << i << " # " << d << " - " << c << std::endl;
+        for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(jacobian,i); it; ++it)
         {
-            for (SparseMatrix::InnerIterator it(m_jacobian,k); it; ++it)
-            {
-                it.valueRef() *= d;
-                if (it.row() == it.col()) it.valueRef() += c;
-            }
+            it.valueRef() *= d;
+            if (it.row() == it.col()) it.valueRef() += c;
         }
     }
+    //std::cout << "Reformulated Jacobian " << std::endl;
+    //std::cout << jacobian.toDense() << std::endl;
 }
 
 // ================================= //
 //                                   //
 //     Solving the linear system     //
 //     -------------------------     //
 //                                   //
 // ================================= //
 
 template <class program_t>
-MiCPParabolicSolverReturnCode MiCPParabolicSolver<program_t>::compute_search_direction(Eigen::VectorXd& update)
+MiCPParabolicSolverReturnCode MiCPParabolicSolver<program_t>::compute_search_direction(Eigen::VectorXd& update,
+                                                                                       ParabolicVariables& x,
+                                                                                       double dt)
 {
-    Eigen::SparseLU<Eigen::SparseMatrix<double, Eigen::RowMajor>,
-                    Eigen::NaturalOrdering<int>> solver;
-    solver.compute(m_jacobian);
+    Eigen::SparseMatrix<double, Eigen::RowMajor> jacobian;
+    setup_jacobian(x, dt, jacobian);
+    Eigen::SparseMatrix<double, Eigen::ColMajor> new_jacobian(jacobian);
+
+    //Eigen::saveMarket(jacobian, "row_major_matrix.mtx");
+    //Eigen::saveMarketVector(-m_phi_residuals, "right_hand_side.mtx");
+    //Eigen::SparseQR<Eigen::SparseMatrix<double, Eigen::RowMajor>,
+    //      Eigen::COLAMDOrdering<int>> solver(jacobian);
+    //std::cout << new_jacobian.toDense() << std::endl;
+    //std::cout << m_phi_residuals << std::endl;
+    Eigen::SparseQR<Eigen::SparseMatrix<double, Eigen::ColMajor>,
+           Eigen::COLAMDOrdering<int>> solver(new_jacobian);
+    //Eigen::BiCGSTAB<Eigen::SparseMatrix<double, Eigen::ColMajor>> solver(new_jacobian);
+    //solver.compute(jacobian);
     m_gradient_step_taken = false;
     if (solver.info() != Eigen::Success)
     {
-        WARNING << "Failed to factorize the jacobian !";
+        ERROR << "Failed to factorize the jacobian !";
         return MiCPParabolicSolverReturnCode::FailedDecomposition;
     }
-    update = solver.solve(- m_phi_residuals);
-    m_gradient = m_jacobian.transpose()*m_phi_residuals;
-    const double condition = - (get_options().factor_descent_direction *
-                                std::pow(update.norm(), get_options().power_descent_direction));
-    if (m_gradient.dot(update) > condition)
+    //solver.setMaxIterations(100);
+    update = solver.solve(-m_phi_residuals);
+    if (solver.info() != Eigen::Success)
     {
-        WARNING << "Gradient step taken !";
-        get_performance().gradient_step_taken = true;
-        update = - m_gradient;
+        ERROR << "Failed to solve the system !";
+        return MiCPParabolicSolverReturnCode::FailedDecomposition;
     }
+    WARNING << "Linear system Solution : " << (new_jacobian*update + m_phi_residuals).norm();
+    //std::cout << solver.iterations() << " - " << solver.error() << std::endl;
+    m_gradient = new_jacobian.transpose()*m_phi_residuals;
+    WARNING << "Update norm : " << update.norm();
+    //std::cout << "Update :" << std::endl << update << std::endl;
+    const double condition = - (get_options().factor_descent_direction *
+                                std::pow(update.norm(), get_options().power_descent_direction));
+    WARNING << "Gradient step condition : " << m_gradient.dot(update) << " compared to : " << condition;
+//    if (m_gradient.dot(update) > condition)
+//    {
+//        WARNING << "Gradient step taken !";
+//        get_performance().gradient_step_taken = true;
+//        update = - m_gradient;
+//    }
+    //std::cout << update << std::endl;
     return MiCPParabolicSolverReturnCode::NotConvergedYet;
 }
 
 template <class program_t>
 double MiCPParabolicSolver<program_t>::is_step_too_long(Eigen::VectorXd& update)
 {
     double steplength = update.norm();
+    //std::cout << "update : " << steplength << std::endl;
     if (steplength > get_options().maximum_step_length)
     {
         get_performance().maximum_step_taken = true;
         update  = get_options().maximum_step_length / steplength * update;
         steplength = get_options().maximum_step_length;
     }
+    //0std::cout << "update : " << steplength << std::endl;
     return steplength;
 }
 
 // ================== //
 //                    //
 //     Linesearch     //
 //     ----------     //
 //                    //
 // ================== //
 
 template <class program_t>
 int MiCPParabolicSolver<program_t>::linesearch(Eigen::VectorXd& update,
                                                ParabolicVariables& x,
                                                const Eigen::VectorXd &predictor,
                                                double dt)
 {
     // References
     // ----------
     //    - Algo A6.3.1 : Dennis and Schnabel (1983)
     //    - Munson et al. (2001)
     //    - Facchinei (2003)
     //    - Nocedal & Wrigth (2006)
     DEBUG << "Linesearch";
     ParabolicVariables xp(x); // copy
     Eigen::VectorXd new_res(get_neq());
     double fcp;
     WARNING << "Norm update : " << update.norm();
     get_performance().maximum_step_taken = false;
     int retcode = 2;
     const double alpha = 1e-4;
     double newtlen = is_step_too_long(update);
     double init_slope = m_gradient.dot(update);
     double rellength = update.cwiseAbs().maxCoeff();
     double minlambda = get_options().step_tolerance / rellength;
-    double lambda = get_program()->maximum_lambda(x, update);
-    double lambda_prev = lambda;
     double alphadt = get_options().alpha*dt;
+    double lambda = get_program()->maximum_lambda(x, update, alphadt);
+    double lambda_prev = lambda;
+
+    //std::cout << "Max lambda : " << lambda << std::endl;
 
     // non monotone linesearch
     // =======================
 
     double merit_value = 0.5*m_phi_residuals.squaredNorm();
     // new residual
+
+//    std::cout << "Before " << std::endl;
+//    int ndf = m_program->get_ndf();
+//    for (int i=0; i<m_program->get_ndf(); ++i)
+//    {
+//    for (int node=0; node<5; ++node)
+//    {
+//        std::cout << xp.displacement(i+ndf*node) << "\t";
+//    }
+//    std::cout << std::endl;
+//    }
+
+
     get_program()->update_variables(xp, predictor, update, lambda, alphadt);
+
+//    std::cout << "After" << std::endl;
+//    for (int i=0; i<m_program->get_ndf(); ++i)
+//    {
+//    for (int node=0; node<5; ++node)
+//    {
+//        std::cout << xp.displacement(i+ndf*node) << "\t";
+//    }
+//    std::cout << std::endl;
+//    }
+
     compute_residuals(xp, new_res);
     reformulate_residuals_inplace(xp, new_res);
+
+    //std::cout << "xp : " << xp.displacement << std::endl;
+    //std::cout << "xp : " << xp.velocity << std::endl;
+    //std::cout << "new_res : " << new_res << std::endl;
+
     fcp = 0.5*new_res.squaredNorm();
 
     // Skip linesearch if enough progress is done
     // ------------------------------------------
     if (fcp < get_options().coeff_accept_newton_step *merit_value)
     {
+        WARNING << "Skip linesearch ";
         if (m_max_merits.size() > 0) m_max_merits[m_max_merits.size()-1] = merit_value;
         else m_max_merits.push_back(merit_value);
         x = xp;
         return 0;
     }
     // Select the merit value of reference
     // -----------------------------------
-    double mmax = merit_value;
-    if (m_max_merits.size() > 0)
-    {
-        mmax = m_max_merits[m_max_merits.size()-1];
-        // check for cycling
-        // - - - - - - - - -
-        if (    m_max_merits.size() ==4
-            &&  std::fabs(mmax - merit_value) < 1e-3*get_options().residuals_tolerance)
-        {
-            WARNING << "Cycling has been detected by the linesearch - Taking the full Newton step";
-            x = xp;
-            update = lambda*update;
-            return 3;
-        }
-    }
-    if (m_max_merits.size() < 4)
-    {
-        m_max_merits.push_back(merit_value);
-        if (merit_value < mmax) merit_value = (3*merit_value + mmax)/4;
-    }
-    else if (merit_value < mmax)
-    {
-        m_max_merits[3] = merit_value;
-        merit_value = mmax;
-    }
-    if (get_performance().gradient_step_taken)
-    {
-        merit_value *= 100;
-    }
+//    double mmax = merit_value;
+//    if (m_max_merits.size() > 0)
+//    {
+//        mmax = m_max_merits[m_max_merits.size()-1];
+//        // check for cycling
+//        // - - - - - - - - -
+//        if (    m_max_merits.size() ==4
+//            &&  std::fabs(mmax - merit_value) < 1e-3*get_options().residuals_tolerance)
+//        {
+//            WARNING << "Cycling has been detected by the linesearch - Taking the full Newton step";
+//            x = xp;
+//            update = lambda*update;
+//            return 3;
+//        }
+//    }
+//    if (m_max_merits.size() < 4)
+//    {
+//        m_max_merits.push_back(merit_value);
+//        if (merit_value < mmax) merit_value = (3*merit_value + mmax)/4;
+//    }
+//    else if (merit_value < mmax)
+//    {
+//        m_max_merits[3] = merit_value;
+//        merit_value = mmax;
+//    }
+    //if (get_performance().gradient_step_taken)
+    //{
+    //    merit_value *= 100;
+    //}
     // The linesearch
     // --------------
     double fc = merit_value;
     double fcp_prev;
     int cnt = 0;
     do
     {
         fcp = 0.5*new_res.squaredNorm();
-        if (fcp <= fc - std::min(-alpha*lambda*init_slope,(1-alpha)*fc)) //pg760 Fachinei2003
+        //WARNING << "cnt : " <<cnt << " - lambda : " << lambda << " # " << fc << " # " << fcp << " # " << init_slope;
+        if (fcp <= fc + alpha*lambda*init_slope)//- std::min(-alpha*lambda*init_slope,(1-alpha)*fc)) //pg760 Fachinei2003
         {
             retcode = 0;
             if (lambda ==1 and (newtlen > 0.99 * get_options().maximum_step_length)) {
                 get_performance().maximum_step_taken = true;
             }
             break;
         }
         else if (lambda < minlambda)
         {
             retcode = 1;
             break;
         }
         else
         {
+            //WARNING << "denom : " << fcp - fc -init_slope;
             // Select a new step length
             // - - - - - - - - - - - -
             double lambdatmp;
             if (cnt == 0) { // only a quadratic at the first
                 lambdatmp = - init_slope / (2*(fcp - fc -init_slope));
             }
             else
             {
                 const double factor = 1 /(lambda - lambda_prev);
                 const double x1 = fcp - fc - lambda*init_slope;
                 const double x2 = fcp_prev - fc - lambda_prev*init_slope;
 
                 const double a = factor * ( x1/(lambda*lambda) - x2/(lambda_prev*lambda_prev));
                 const double b = factor * ( -x1*lambda_prev/(lambda*lambda) + x2*lambda/(lambda_prev*lambda_prev));
 
                 if (a == 0)
                 { // cubic interpolation is in fact a quadratic interpolation
                     lambdatmp = - init_slope/(2*b);
                 }
                 else
                 {
                     const double disc = b*b-3*a*init_slope;
                     lambdatmp = (-b+std::sqrt(disc))/(3*a);
                 }
                 if (lambdatmp > 0.5*lambda ) lambdatmp = 0.5*lambda;
             }
+            //WARNING << "lambdatmp : " << lambdatmp;
             lambda_prev = lambda;
             fcp_prev = fcp;
             if (lambdatmp < 0.1*lambda) {
                 lambda = 0.1 * lambda;
             } else {
                 lambda = lambdatmp;
             }
         }
         get_program()->update_variables(xp, predictor, update, lambda, alphadt);
         compute_residuals(xp, new_res);
         reformulate_residuals_inplace(xp, new_res);
         ++cnt;
     } while(retcode == 2 and cnt < 100);
     WARNING << "Lambda : " << lambda;
     if (cnt == 100)
     {
         ERROR << "Too much linesearch iterations ! We stop";
     }
     get_program()->update_variables(x, predictor, update, lambda, alphadt);
     return retcode;
 }
 
 } // end namespace micpsolver
 } // end namespace reactmicp
 } // end namespace specmicp
diff --git a/src/reactmicp/micpsolvers/micpsolver_structs.hpp b/src/reactmicp/micpsolvers/micpsolver_structs.hpp
index a457908..d21499d 100644
--- a/src/reactmicp/micpsolvers/micpsolver_structs.hpp
+++ b/src/reactmicp/micpsolvers/micpsolver_structs.hpp
@@ -1,82 +1,82 @@
 #ifndef SPECMICP_REACTMICP_MICPSOLVER_MICPSOLVERSTRUCTS_HPP
 #define SPECMICP_REACTMICP_MICPSOLVER_MICPSOLVERSTRUCTS_HPP
 
 namespace specmicp {
 namespace reactmicp {
 namespace micpsolver  {
 
 //! options for the ReactMiCP solver
 struct MiCPParabolicSolverOptions
 {
     double alpha;
 
     double penalization_factor;
     double factor_descent_direction;
     double power_descent_direction;
 
     double step_tolerance;
     double residuals_tolerance;
     double threshold_stationary_point;
 
     double maximum_step_length;
     double coeff_accept_newton_step;
     double maximum_iterations;
     double maximum_iterations_at_maximum_step;
 
     MiCPParabolicSolverOptions():
         alpha(1.0),
         penalization_factor(0.8),
-        factor_descent_direction(0.5),
-        power_descent_direction(2.0),
-        step_tolerance(1e-6),
+        factor_descent_direction(1.5),
+        power_descent_direction(1.0),
+        step_tolerance(1e-8),
         residuals_tolerance(1e-6),
         threshold_stationary_point(1e-4),
-        maximum_step_length(1e3),
+        maximum_step_length(100.0),
         coeff_accept_newton_step(0.9),
-        maximum_iterations(200),
+        maximum_iterations(50),
         maximum_iterations_at_maximum_step(50)
     {}
 
 };
 
 enum class MiCPParabolicSolverReturnCode
 {
     UnknownError = -10,
     FailedDecomposition = -5,
     MaxStepTakenTooManyTimes = -3,
     MaxIterations = -2,
     StationaryPoint = -1,
     NotConvergedYet = 0,
     ResidualMinimized = 1,
     ErrorMinimized = 2
 };
 
 struct MiCPParabolicSolverPerformance
 {
     bool gradient_step_taken;
     bool maximum_step_taken;
     int nb_iterations;
     int nb_maximum_step_taken;
     int nb_consecutive_maximum_step_taken;
     int nb_call_residuals;
     int nb_call_jacobian;
     MiCPParabolicSolverReturnCode return_code;
 
     MiCPParabolicSolverPerformance():
         gradient_step_taken(false),
         maximum_step_taken(false),
         nb_iterations(0),
         nb_maximum_step_taken(0),
         nb_consecutive_maximum_step_taken(0),
         nb_call_residuals(0),
         nb_call_jacobian(0),
         return_code(MiCPParabolicSolverReturnCode::NotConvergedYet)
     {}
 };
 
 } // end namespace micpsolver
 } // end namespace reactmicp
 } // end namespace specmicp
 
 
 #endif // SPECMICP_REACTMICP_MICPSOLVER_MICPSOLVERSTRUCTS_HPP
diff --git a/src/reactmicp/micpsolvers/parabolicprogram.hpp b/src/reactmicp/micpsolvers/parabolicprogram.hpp
index 8053aed..e656358 100644
--- a/src/reactmicp/micpsolvers/parabolicprogram.hpp
+++ b/src/reactmicp/micpsolvers/parabolicprogram.hpp
@@ -1,113 +1,113 @@
 /*-------------------------------------------------------
 
  - Module : reactmicp
  - File : common.hpp
  - Author : Fabien Georget
 
 Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
     * Neither the name of the Princeton University nor the
       names of its contributors may be used to endorse or promote products
       derived from this software without specific prior written permission.
 
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 ---------------------------------------------------------*/
 
 #ifndef SPECMICP_REACTMICP_PARABOLICPROGRAM_HPP
 #define SPECMICP_REACTMICP_PARABOLICPROGRAM_HPP
 
 #include <reactmicp/common.hpp>
 
 namespace specmicp {
 namespace reactmicp {
 
 namespace solvers {
 
 
 //! \brief A abstract base class for a parabolic program
 //!
 template <class Derived>
 class ParabolicProgram
 {
 public:
     //! \brief Return the number of degree of freedoms
     ind_t get_neq()
     {
         return derived()->get_neq();
     }
 
 
     //! \brief Return the derived class
     Derived& derived() {return static_cast<Derived*>(*this);}
 
     //! \brief Return the residuals
     void get_residuals(const ParabolicVariables& variable,
                        Vector& residual)
     {
         return derived()->get_residuals(variable, residual);
     }
 
     //! brief Return the jacobian
     void get_jacobian(const ParabolicVariables& variable,
                       Vector& residual,
                       SparseMatrix& jacobian)
     {
         return derived()->get_jacobian(variable, residual, jacobian);
     }
 
     //! \brief Return the index of the MiCP equations
     std::vector<bool>& micp_index()
     {
         return derived()->micp_index();
     }
 
     double maximum_lambda(const Eigen::VectorXd& x,
                           const Eigen::VectorXd& update);
 
     //! \brief Called at the beginning of an iteration
     //!
     //! @return A boolean indicating if the system can have converged
     //!
     //! Return true by default
     bool hook_start_iteration(const ParabolicVariables& x, double norm_residual) { return true;}
 
     //! \brief Updates the variables
     void update_variables(ParabolicVariables& x,
                           const Eigen::VectorXd& predictor,
                           const Eigen::VectorXd& update,
                           double factor
                           );
 
     //! \brief Set the predictor
-    void set_predictor(const ParabolicVariables& x,
+    void set_predictor(ParabolicVariables& x,
                        Eigen::VectorXd& predictor,
                        double alpha,
                        double dt);
 
     //! \brief Return the maximum lambda allowed for the linesearch
-    double maximum_lambda(const ParabolicVariables& x, const Eigen::VectorXd& update) {return 1;}
+    double maximum_lambda(const ParabolicVariables& x, const Eigen::VectorXd& update, double alphadt) {return 1.0;}
 };
 
 } // end namespace solvers
 } // end namespace reactmicp
 } // end namespace specmicp
 
 #endif // SPECMICP_REACTMICP_PARABOLICPROGRAM_HPP
diff --git a/src/reactmicp/systems/diffusion/diffusion.cpp b/src/reactmicp/systems/diffusion/diffusion.cpp
index 937f653..39805b3 100644
--- a/src/reactmicp/systems/diffusion/diffusion.cpp
+++ b/src/reactmicp/systems/diffusion/diffusion.cpp
@@ -1,373 +1,521 @@
 /*-------------------------------------------------------
 
  - Module : reactmicp/systems/diffusion
  - File : diffusion.cpp
  - Author : Fabien Georget
 
 Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
     * Neither the name of the Princeton University nor the
       names of its contributors may be used to endorse or promote products
       derived from this software without specific prior written permission.
 
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 ---------------------------------------------------------*/
 
 #include "diffusion.hpp"
 #include "physics/laws.hpp"
 #include "utils/log.hpp"
 #include <iostream>
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 
 
 DiffusionProgram::DiffusionProgram(std::shared_ptr<mesh::Mesh1D> themesh,
                  std::shared_ptr<database::DataContainer> thedatabase,
                  std::shared_ptr<DiffusionParameter> parameters,
                  diffusion::ListBoundaryCondition& list_bcs,
-                 EquilibriumState& initialstate,
+                 const EquilibriumState& initialstate,
                  Variables& var):
     DiffusionProgram(themesh, thedatabase, parameters, list_bcs, var)
 {
-    for (ind_t node: themesh->range_nodes())
-    {
-        if (not m_ideq.node_has_bc(node))
-        {
-            for (ind_t component: thedatabase->range_aqueous_component())
-            {
-                m_ideq.component_concentration(node, component, var) = std::log10(initialstate.molality_component(component));
-                m_ideq.mobile_total_concentration(node, component, var) = initialstate.total_aqueous_concentration_component(component);
-                m_second.loggamma_component(node, component) = std::log10(initialstate.activity_coefficient_component(component));
-            }
-            for (ind_t aqueous: thedatabase->range_aqueous())
-            {
-                m_second.secondary_concentration(node, aqueous) = initialstate.molality_secondary()[aqueous];
-                m_second.loggamma_secondary(node, aqueous) = std::log10(initialstate.activity_coefficient_secondary(aqueous));
-            }
-            for (ind_t mineral: thedatabase->range_mineral())
-            {
-                m_ideq.mole_mineral(node, mineral, var) = initialstate.moles_mineral(mineral);
-            }
-        }
-    }
+    initialize_no_bc_with_init(initialstate, var);
 }
 
 // ================================== //
 //                                    //
 //             Residuals              //
 //                                    //
 // ================================== //
 
 void DiffusionProgram::get_residuals(const Variables& variable,
                                      Vector& residual)
 {
-    m_is_micp = std::vector<bool>(m_ideq.get_neq(), false);
+    m_is_micp = std::vector<bool>(get_neq(), false);
+    residual = Eigen::VectorXd::Zero(get_neq());
     residual_transport(variable, residual);
     residual_speciation(variable, residual);
 }
 
 // Transport
 // =========
 
 void DiffusionProgram::residual_transport(const Variables& variable,
                                           Vector& residual)
 {
     for (ind_t element: m_mesh->range_elements())
     {
         element_residual_transport(element, variable, residual);
     }
 }
 
+void DiffusionProgram::element_residual_transport_component(
+        ind_t element,
+        int component,
+        const Variables& variable,
+        Vector& element_residual)
+{
+
+    Eigen::Matrix2d mass, jacob;
+    Eigen::Vector2d velocity, displacement;
+    double mass_coeff = -(m_param->density_water()*
+                          m_param->porosity(0)*
+                          m_mesh->crosssection()*m_mesh->get_dx(element)/2);
+    mass << 1, 0, 0, 1;
+    mass *= mass_coeff;
+    double flux_coeff = -( m_mesh->crosssection() / m_mesh->get_dx(element)
+                           * m_param->porosity(0)
+                           * m_param->diffusion_coefficient(element)
+                           *m_param->density_water()
+                           );
+    jacob << 1, -1, -1, 1;
+    jacob *= flux_coeff;
+    velocity << variable.velocity(m_ideq.get_dof_diffusion(m_mesh->get_node(element,0), component)),
+            variable.velocity(m_ideq.get_dof_diffusion(m_mesh->get_node(element,1), component));
+    displacement << m_ideq.mobile_total_concentration(m_mesh->get_node(element,0), component, variable),
+            m_ideq.mobile_total_concentration(m_mesh->get_node(element,1), component, variable);
+
+    element_residual = mass*velocity+jacob*displacement;
+
+    for (int en: m_mesh->range_nen())
+    {
+        element_residual(en) += nodal_mineral_transient_term_transport(
+                    m_mesh->get_node(element, en), component, variable)/2;
+    }
+}
+
 void DiffusionProgram::element_residual_transport(ind_t element,
                                                   const Variables& variable,
                                                   Vector& residual)
 {
+    Eigen::VectorXd element_residual(2);
     for (ind_t component: m_data->range_aqueous_component())
     {
-        double face_flux = ( m_mesh->crosssection()/m_mesh->get_dx(element)
-                             * m_param->diffusion_coefficient(element)*
-                             (
-                                 m_ideq.mobile_total_concentration(m_mesh->get_node(element,0), component, variable)
-                               - m_ideq.mobile_total_concentration(m_mesh->get_node(element,1), component, variable)
-                              )
-                            );
-        for (ind_t indnode: m_mesh->range_nen())
+        element_residual_transport_component(element, component, variable, element_residual);
+
+        for (int en: m_mesh->range_nen())
         {
-            ind_t node = m_mesh->get_node(element, indnode);
-            ind_t id = m_ideq.id_equation_diffusion(node, component);
-            if (id == no_equation) continue;
-            // mineral
-            double res = nodal_mineral_transient_term_transport(node, component, variable);
-            // transient term
-            res -= m_param->porosity(node)*variable.velocity(m_ideq.get_dof_diffusion(node, component));
-            res *= m_mesh->crosssection()*m_mesh->get_dx(element)/2;
-            // diffusion
-            res += ((node==0)?1:-1)*face_flux;
-            // update
-            residual(id) += res;
+            const ind_t node = m_mesh->get_node(element, en);
+            const ind_t id = m_ideq.id_equation_diffusion(node, component);
+            if (id != no_equation) {residual.coeffRef(id) += element_residual(en);}
         }
     }
 }
 
 double DiffusionProgram::nodal_mineral_transient_term_transport(ind_t node,
                                                                 ind_t component,
                                                               const Variables& variable)
 {
     double mineral_term = 0;
     for (ind_t mineral: m_data->range_mineral())
     {
         if (m_data->nu_mineral(mineral, component) == 0 ) continue;
         mineral_term -= ( m_data->nu_mineral(mineral, component) *
                           variable.velocity(m_ideq.get_dof_mineral(node, mineral))
                          );
     }
-    mineral_term *= 1.0/m_mesh->cell_volume(node);
+    //mineral_term *= 1.0/2.0; //(2.0/m_mesh->cell_volume(node));
     return mineral_term;
 }
 
 // Speciation
 // ==========
 
 void DiffusionProgram::residual_speciation(const Variables& variable,
                                            Vector& residual)
 {
     for (ind_t node: m_mesh->range_nodes())
     {
         nodal_residual_massbalance(node, variable, residual);
         nodal_residual_mineral(node, variable, residual);
     }
 }
 
 void DiffusionProgram::nodal_residual_massbalance(ind_t node,
                                                    const Variables& variable,
                                                    Vector& residual)
 {
     for (ind_t component: m_data->range_aqueous_component())
     {
         ind_t id = m_ideq.id_equation_massbalance(node, component);
         if (id != no_equation)
         {
             double sum = pow10(m_ideq.component_concentration(node, component, variable));
             for (ind_t aqueous: m_data->range_aqueous())
             {
                 if (m_data->nu_aqueous(aqueous, component) != 0)
                 {
                     sum += m_data->nu_aqueous(aqueous, component)*m_second.secondary_concentration(node, aqueous);
                 }
             }
-            residual(id) += m_ideq.mobile_total_concentration(node, component, variable)/m_param->density_water() - sum ;
+            residual(id) += sum - m_ideq.mobile_total_concentration(node, component, variable);///m_param->density_water();
         }
     }
 }
 
 void DiffusionProgram::nodal_residual_mineral(ind_t node,
                                               const Variables& variable,
                                               Vector& residual)
 {
     for (ind_t mineral: m_data->range_mineral())
     {
         ind_t id  = m_ideq.id_equation_mineral(node, mineral);
         if (id != no_equation)
         {
             double sum = m_data->logk_mineral(mineral);
             for (ind_t component: m_data->range_aqueous_component())
             {
                 sum -= m_data->nu_mineral(mineral, component) * (
                             m_ideq.component_concentration(node, component, variable)
                             + m_second.loggamma_component(node, component));
             }
             residual(id) = sum;
             m_is_micp[id] = true;
         }
     }
 }
 
 // ================================== //
 //                                    //
 //            Jacobian                //
 //                                    //
 // ================================== //
 
-void DiffusionProgram::get_jacobian(const Variables &variable,
+void DiffusionProgram::get_jacobian(Variables &variable,
                                     Vector &residual,
                                     SparseMatrix &jacobian,
                                     double alphadt)
 {
     list_triplet_t jacob;
     const ind_t ncomp = m_data->nb_component-1;
     const ind_t nmin = m_data->nb_mineral;
     const ind_t estimation = m_mesh->nnodes()*(ncomp*(3+ncomp) + ncomp*(ncomp+nmin)+nmin+ncomp*2);
     jacob.reserve(estimation);
     jacobian_transport(variable, residual, jacob, alphadt);
     jacobian_speciation(variable, jacob, alphadt);
     jacobian = SparseMatrix(get_neq(), get_neq());
     jacobian.setFromTriplets(jacob.begin(), jacob.end());
 }
 
 // Transport
 // =========
 
-void DiffusionProgram::jacobian_transport(const Variables& variable,
+void DiffusionProgram::jacobian_transport(Variables& variable,
                                           Vector& residual,
                                           list_triplet_t& jacobian,
                                           double alphadt)
 {
     for (ind_t element: m_mesh->range_elements())
     {
-        element_jacobian_transport(element, variable, residual, jacobian, alphadt);
+        element_jacobian_transport_fd(element, variable, residual, jacobian, alphadt);
+    }
+}
+
+
+void DiffusionProgram::element_jacobian_transport_fd(ind_t element,
+                                                  Variables& variable,
+                                                  Vector& residual,
+                                                  list_triplet_t& jacobian,
+                                                  double alphadt)
+{
+
+    const double eps = 1e-6;
+
+    for (int component: m_data->range_aqueous_component())
+    {
+        Eigen::VectorXd element_residual_orig(Eigen::VectorXd::Zero(2));
+        element_residual_transport_component(element, component, variable, element_residual_orig);
+
+        //ind_t idl = m_ideq.id_equation_diffusion(m_mesh->get_node(element, 0), component);
+        //ind_t idr = m_ideq.id_equation_diffusion(m_mesh->get_node(element, 1), component);
+
+        for (int en: m_mesh->range_nen())
+        {
+            Eigen::VectorXd element_residual(Eigen::VectorXd::Zero(2));
+
+            const ind_t node = m_mesh->get_node(element, en);
+            const ind_t idc = m_ideq.id_equation_diffusion(node, component);
+            const ind_t dof = m_ideq.get_dof_diffusion(node, component);
+
+            if (idc == no_equation) continue;
+
+            const double tmp_v = variable.velocity(dof);
+            const double tmp_d = variable.displacement(dof);
+
+            double h = eps*std::abs(tmp_v);
+            if (h == 0) h = eps;
+            variable.velocity(dof) = tmp_v + h;
+            h = variable.velocity(dof) - tmp_v;
+
+            variable.displacement(dof) = tmp_d + alphadt*h;
+
+            element_residual_transport_component(element, component, variable, element_residual);
+            variable.velocity(dof) = tmp_v;
+            variable.displacement(dof) = tmp_d;
+
+            for (int enr: m_mesh->range_nen())
+            {
+                const ind_t noder = m_mesh->get_node(element, enr);
+                const ind_t idr = m_ideq.id_equation_diffusion(noder, component);
+
+                if (idr == no_equation) continue;
+                jacobian.push_back(triplet_t(idr, idc, (element_residual(enr) - element_residual_orig(enr))/h));
+            }
+            // mineral -> not using finite difference
+            for (ind_t mineral: m_data->range_mineral())
+            {
+                const ind_t idm = m_ideq.id_equation_mineral(node, mineral);
+                if (idm != no_equation)
+                    jacobian.push_back(triplet_t(idc, idm, -m_data->nu_mineral(mineral, component)/2.0));
+            }
+        }
+
     }
 }
 
 void DiffusionProgram::element_jacobian_transport(ind_t element,
                                                   const Variables& variable,
                                                   Vector& residual,
                                                   list_triplet_t& jacobian,
                                                   double alphadt)
 {
-    double deriv = alphadt*( m_mesh->crosssection()/m_mesh->get_dx(element)
-                     * m_param->diffusion_coefficient(element));
+    double deriv = -alphadt*( m_mesh->crosssection() /m_mesh->get_dx(element)
+                     * m_param->diffusion_coefficient(element)*m_param->porosity(0)
+                     * m_param->density_water());
     for (ind_t component: m_data->range_aqueous_component())
     {
         ind_t idl = m_ideq.id_equation_diffusion(m_mesh->get_node(element, 0), component);
         ind_t idr = m_ideq.id_equation_diffusion(m_mesh->get_node(element, 1), component);
         if (idl != no_equation)
         {
             jacobian.push_back(triplet_t(idl, idl, deriv));
             if (idr != no_equation) jacobian.push_back(triplet_t(idl, idr, -deriv));
         }
         if (idr != no_equation)
         {
             jacobian.push_back(triplet_t(idr, idr, deriv));
             if (idl != no_equation) jacobian.push_back(triplet_t(idr, idl, -deriv));
         }
         // Mass Matrix
-        double mass_coeff = (m_param->porosity(m_mesh->get_node(element, 0)) + m_mesh->get_node(element, 1))/2;
-        mass_coeff *= m_mesh->crosssection()*m_mesh->get_dx(element)/2;
+        double mass_coeff = -(m_param->porosity(m_mesh->get_node(element, 0)));
+        mass_coeff *= m_param->density_water()*m_mesh->crosssection()*m_mesh->get_dx(element)/2;
         if (idl != no_equation) jacobian.push_back(triplet_t(idl, idl, mass_coeff));
         if (idr != no_equation) jacobian.push_back(triplet_t(idr, idr, mass_coeff));
+        // mineral
+//        for (ind_t mineral: m_data->range_mineral())
+//        {
+//            ind_t idm_0 = m_ideq.id_equation_mineral(m_mesh->get_node(element, 0), mineral);
+//            ind_t idm_1 = m_ideq.id_equation_mineral(m_mesh->get_node(element, 1), mineral);
+//            if (idm_0 != no_equation and idl != no_equation)
+//                jacobian.push_back(triplet_t(idl, idm_0,
+//                             -m_data->nu_mineral(mineral, component)/2.0));
+//            if (idm_1 != no_equation and idr != no_equation)
+//                jacobian.push_back(triplet_t(idr, idm_1,
+//                             -m_data->nu_mineral(mineral, component)/2.0));
+//        }
     }
 }
 
 // Speciation
 // ==========
 
 void DiffusionProgram::jacobian_speciation(const Variables& variable,
                                            list_triplet_t& jacobian,
                                            double alphadt)
 {
     for (ind_t node: m_mesh->range_nodes())
     {
         nodal_jacobian_massbalance(node, variable, jacobian, alphadt);
         nodal_jacobian_mineral(node, variable, jacobian, alphadt);
 
     }
 }
 
 void DiffusionProgram::nodal_jacobian_massbalance(ind_t node,
                                                  const Variables& variable,
                                                  list_triplet_t& jacobian,
                                                   double alphadt)
 {
     const double logten = std::log(10.0);
     for (ind_t component: m_data->range_aqueous_component())
     {
         const ind_t idp = m_ideq.id_equation_massbalance(node, component);
         if (idp == no_equation) continue;
+        // total concentration
+        const ind_t idc = m_ideq.id_equation_diffusion(node, component);
+        if (idc != no_equation)
+        {
+            jacobian.push_back(triplet_t(idp, idc, -alphadt)); // /m_param->density_water()));
+        }
         for (ind_t k: m_data->range_aqueous_component())
         {
-            // total concentration
-            const ind_t idc = m_ideq.id_equation_diffusion(node, component);
-            if (idc != no_equation)
-            {
-                jacobian.push_back(triplet_t(idp, idc, alphadt/m_param->density_water()));
-            }
             // concentration
             const ind_t ids = m_ideq.id_equation_massbalance(node, component);
             if (ids == no_equation) continue;
             double tmp_iip = 0;
-            if (k == component) tmp_iip -= logten*pow10(m_ideq.component_concentration(node, component, variable));
+            if (k == component) tmp_iip += logten*pow10(m_ideq.component_concentration(node, component, variable));
             for (ind_t j: m_data->range_aqueous())
             {
-                tmp_iip -= logten*m_data->nu_aqueous(j, component)*m_data->nu_aqueous(j, k)*m_second.secondary_concentration(node, j);
+                tmp_iip += (logten*m_data->nu_aqueous(j, component)*
+                            m_data->nu_aqueous(j, k)*m_second.secondary_concentration(node, j));
             }
             jacobian.push_back(triplet_t(idp, ids, alphadt*tmp_iip));
         }
     }
 }
 
 void DiffusionProgram::nodal_jacobian_mineral(ind_t node,
                                                  const Variables& variable,
                                                  list_triplet_t& jacobian,
                                               double alphadt)
 {
     for (ind_t mineral: m_data->range_mineral())
     {
         const ind_t idm = m_ideq.id_equation_mineral(node, mineral);
         if (idm == no_equation) continue;
         for (ind_t component: m_data->range_aqueous_component())
         {
             const ind_t idc = m_ideq.id_equation_massbalance(node, component);
             if (idc == no_equation) continue;
             jacobian.push_back(triplet_t(idm, idc, -alphadt*m_data->nu_mineral(mineral, component)));
         }
         // make place for an element
         jacobian.push_back(triplet_t(idm, idm, 0));
     }
 }
 
 // Variables
 // =========
 
 void DiffusionProgram::update_variables(ParabolicVariables& x,
                       const Eigen::VectorXd& predictor,
                       const Eigen::VectorXd& update,
                       double factor,
                       double alpha_dt
                       )
 {
     for (ind_t node: m_mesh->range_nodes())
     {
         for (ind_t edof=0; edof<get_ndf(); ++edof)
         {
             if (m_ideq.id_equation(node, edof) == no_equation) continue;
             const ind_t dof = m_ideq.get_dof(node, edof);
             const ind_t id_eq = m_ideq.id_equation(node, edof);
             x.velocity(dof) += factor*update(id_eq);
             x.displacement(dof) = predictor(dof) + alpha_dt * x.velocity(dof);
         }
     }
 }
 
-void DiffusionProgram::set_predictor(const ParabolicVariables& x,
+void DiffusionProgram::set_predictor(ParabolicVariables& x,
                    Eigen::VectorXd& predictor,
                    double alpha,
                    double dt)
 {
     predictor.resize(m_mesh->nnodes()*get_ndf());
     predictor = x.displacement + (1-alpha)*dt*x.velocity;
+    x.velocity.setZero();
     WARNING << "Predictor size : " << predictor.rows();
 }
 
+// compute secondary conc
+bool DiffusionProgram::hook_start_iteration(const Variables& x, double norm_residual) {
+    return m_second.solve_secondary_variables(x);
+}
+
+
+double DiffusionProgram::maximum_lambda(const ParabolicVariables& x, const Eigen::VectorXd& update, double alphadt)
+{
+    double inv_maximum_lambda = 1;
+    for (ind_t node: m_mesh->range_nodes())
+    {
+        for (int mineral: m_data->range_mineral())
+        {
+            if (m_ideq.id_equation_mineral(node, mineral) == no_equation
+                    or x.displacement(m_ideq.get_dof_mineral(node, mineral)) == 0) continue;
+            inv_maximum_lambda = std::max(inv_maximum_lambda,
+                        -alphadt*update(m_ideq.id_equation_mineral(node, mineral))/
+                               (0.9*x.displacement(m_ideq.get_dof_mineral(node, mineral))));
+        }
+    }
+    return 1.0/inv_maximum_lambda;
+}
+
+void DiffusionProgram::initialize_no_bc_with_init(const EquilibriumState& initialstate,
+                                                  ParabolicVariables& var)
+{
+    for (ind_t node: m_mesh->range_nodes())
+    {
+        if (not m_ideq.node_has_bc(node))
+        {
+            for (ind_t component: m_data->range_aqueous_component())
+            {
+                m_ideq.component_concentration(node, component, var) =
+                        std::log10(initialstate.molality_component(component));
+                m_ideq.mobile_total_concentration(node, component, var) =
+                        initialstate.total_aqueous_concentration_component(component);
+                m_second.loggamma_component(node, component) = initialstate.loggamma_component(component);
+            }
+            for (ind_t aqueous: m_data->range_aqueous())
+            {
+                m_second.secondary_concentration(node, aqueous) = initialstate.molality_secondary()[aqueous];
+                m_second.loggamma_secondary(node, aqueous) = initialstate.loggamma_secondary(aqueous);
+            }
+            for (ind_t mineral: m_data->range_mineral())
+            {
+                m_ideq.mole_mineral(node, mineral, var) = initialstate.moles_mineral(mineral);
+            }
+        }
+
+    }
+}
+
+void DiffusionProgram::copy_cp_variables(const ParabolicVariables& x,
+                                         Eigen::VectorXd &moles_mineral)
+{
+    moles_mineral.resize(get_neq());
+    moles_mineral.setZero();
+    for (ind_t node: m_mesh->range_nodes())
+    {
+        for (int mineral: m_data->range_mineral())
+        {
+            const int id_eq =  m_ideq.id_equation_mineral(node, mineral);
+            if (id_eq == no_equation) continue;
+            moles_mineral(id_eq) = m_ideq.mole_mineral(node, mineral, x);
+        }
+    }
+}
+
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
diff --git a/src/reactmicp/systems/diffusion/diffusion.hpp b/src/reactmicp/systems/diffusion/diffusion.hpp
index d899211..378b1bc 100644
--- a/src/reactmicp/systems/diffusion/diffusion.hpp
+++ b/src/reactmicp/systems/diffusion/diffusion.hpp
@@ -1,232 +1,259 @@
 /*-------------------------------------------------------
 
  - Module : reactmicp/systems/diffusion
  - File : diffusion.hpp
  - Author : Fabien Georget
 
 Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
     * Neither the name of the Princeton University nor the
       names of its contributors may be used to endorse or promote products
       derived from this software without specific prior written permission.
 
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 ---------------------------------------------------------*/
 
 #ifndef SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP
 #define SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP
 
 #include <vector>
 #include <memory>
 
 #include "database/data_container.hpp"
 
 #include "reactmicp/common.hpp"
 #include "reactmicp/micpsolvers/parabolicprogram.hpp"
 #include "reactmicp/meshes/mesh1d.hpp"
 
 #include "diffusion_numbering.hpp"
 #include "diffusion_secondary.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 
 //! \brief Parameters for the diffusion system
 class DiffusionParameter
 {
 public:
 
     DiffusionParameter(double diffusion_coefficient, double porosity):
         m_diff(diffusion_coefficient), m_poro(porosity)
     {}
 
     //! \brief Density of water (kg/m^3)
     double density_water() {return 1e3;}
 
     //! \brief Diffusion coefficient (element by element) (m^2/s)
     double diffusion_coefficient(ind_t element) {return m_diff;}
 
     //! \brief Porosity (at a node (function of composition))
     double porosity(ind_t node) {return m_poro;}
 
 private:
     double m_diff;
     double m_poro;
 };
 
 //! \brief Saturated diffusion in reactive porous media
 class DiffusionProgram: public solvers::ParabolicProgram<DiffusionProgram>
 {
 
 public:
     using Variables = ParabolicVariables;
     const int no_equation = -1;
 
     DiffusionProgram(std::shared_ptr<mesh::Mesh1D> themesh,
                      std::shared_ptr<database::DataContainer> thedatabase,
                      std::shared_ptr<DiffusionParameter> parameters,
                      diffusion::ListBoundaryCondition& list_bcs,
                      Variables& var):
         m_mesh(themesh), m_data(thedatabase), m_param(parameters),
         m_ideq(themesh->nnodes(), thedatabase, list_bcs, var),
         m_second(themesh->nnodes(), thedatabase, m_ideq)
     {
 
     }
 
     DiffusionProgram(std::shared_ptr<mesh::Mesh1D> themesh,
                      std::shared_ptr<database::DataContainer> thedatabase,
                      std::shared_ptr<DiffusionParameter> parameters,
                      diffusion::ListBoundaryCondition& list_bcs,
-                     EquilibriumState& initialstate,
+                     const EquilibriumState &initialstate,
                      Variables& var);
 
 
     //! \brief Return the number of degree of freedoms
     ind_t get_neq() const {return m_ideq.get_neq();}
     ind_t get_ndf() const {return m_ideq.get_ndf();}
 
     //! \brief Return the residuals
-    void get_residuals(const Variables& variable,
+    void get_residuals(const ParabolicVariables& variable,
                        Vector& residual);
 
     //! brief Return the jacobian
-    void get_jacobian(const Variables& variable,
+    void get_jacobian(ParabolicVariables& variable,
                       Vector& residual,
                       SparseMatrix& jacobian,
                       double alphadt);
 
     //! \brief Return the index of the MiCP equations
     const std::vector<bool>& micp_index() const { return m_is_micp;}
 
 
     //! \brief Update the variables
     void update_variables(ParabolicVariables& x,
                           const Eigen::VectorXd& predictor,
                           const Eigen::VectorXd& update,
                           double factor,
                           double alpha_dt
                           );
     //! \brief set the predictor
-    void set_predictor(const ParabolicVariables& x,
+    void set_predictor(ParabolicVariables& x,
                        Eigen::VectorXd& predictor,
                        double alpha,
                        double dt);
+
+    //! \brief compute secondary concentrations
+    bool hook_start_iteration(const ParabolicVariables& x, double norm_residual);
+
+    //! \brief Return the maximum lambda allowed for the linesearch
+    double maximum_lambda(const ParabolicVariables& x, const Eigen::VectorXd& update, double alphadt);
+
+    //! \brief copy mineral amount from variables
+    void copy_cp_variables(const ParabolicVariables& x,
+                      Eigen::VectorXd& moles_mineral);
 private:
 
     // Residuals
     // =========
 
     // transport
     // ---------
 
     //! \brief Compute the residuals for the diffusion equations
     void residual_transport(const Variables& variable,
                                               Vector& residual);
 
     //! \brief Contribution of the diffusion operator to the residuals
     void element_residual_transport(ind_t element,
                                     const Variables& variable,
                                     Vector& residual);
 
     //! \brief Contribution of the minerals to the residuals of the diffusion equations
     double nodal_mineral_transient_term_transport(ind_t node,
                                                   ind_t component,
                                                 const Variables& variable);
 
+    void element_residual_transport_component(
+            ind_t element,
+            int component,
+            const Variables& variable,
+            Vector& element_residual);
+
     // speciation
     // ----------
 
     //! \brief Compute the residual for the speciation problem
     void residual_speciation(const Variables& variable,
                                Vector& residual);
 
     //! \brief Compute the residuals of the aqueous mass balances
     void nodal_residual_massbalance(ind_t node,
                                     const Variables& variable,
                                     Vector& residual);
 
     //! \brief Compute the residuals of the mineral speciation
     void nodal_residual_mineral(ind_t node,
                                 const Variables& variable,
                                 Vector& residual);
 
     // Jacobian
     // ========
 
     // transport
     // ---------
 
     //! \brief Contribution of the diffusion equation to the jacobian
-    void jacobian_transport(const Variables& variable,
+    void jacobian_transport(Variables& variable,
                             Vector& residual,
                             list_triplet_t& jacobian,
                             double alphadt);
 
     //! \brief Element by element assembly procedure
     void element_jacobian_transport(ind_t element,
                                     const Variables& variable,
                                     Vector& residual,
                                     list_triplet_t& jacobian,
                                     double alphadt);
+
+    void element_jacobian_transport_fd(ind_t element,
+                                                      Variables &variable,
+                                                      Vector& residual,
+                                                      list_triplet_t& jacobian,
+                                                      double alphadt);
     // speciation
     // ----------
 
     //! \brief Contribution of the speciation problem to the mass balance
     void jacobian_speciation(const Variables& variable,
                              list_triplet_t& jacobian,
                              double alphadt);
 
     //! \brief Contribution of the mass balance at 'node' to the jacobian
     void nodal_jacobian_massbalance(ind_t node,
                                    const Variables& variable,
                                    list_triplet_t& jacobian,
                                     double alphadt);
 
     //! \brief Contribution of the mineral at equilibrium problem to the jacobian
     void nodal_jacobian_mineral(ind_t node,
                                 const Variables& variable,
                                 list_triplet_t& jacobian,
                                 double alphadt);
 
+    // Initialization
+    // ==============
+
+    void initialize_no_bc_with_init(const EquilibriumState& initialstate,
+                                    ParabolicVariables& var);
 
     // Member
     // =======
 
     std::shared_ptr<mesh::Mesh1D> m_mesh;
     std::shared_ptr<database::DataContainer> m_data;
     std::shared_ptr<DiffusionParameter> m_param;
 
     diffusion::DiffusionNumbering m_ideq;
     diffusion::DiffusionSecondaryVariables m_second;
 
     std::vector<bool> m_is_micp;
 
     //SparseVector m_boundary_condition;
 
 };
 
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
 
 #endif // SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP
diff --git a/src/reactmicp/systems/diffusion/diffusion_numbering.cpp b/src/reactmicp/systems/diffusion/diffusion_numbering.cpp
index 4dcce6d..07cf9a0 100644
--- a/src/reactmicp/systems/diffusion/diffusion_numbering.cpp
+++ b/src/reactmicp/systems/diffusion/diffusion_numbering.cpp
@@ -1,151 +1,168 @@
 /*-------------------------------------------------------
 
  - Module : reactmicp/systems/diffusion
  - File : diffusion_numbering.cpp
  - Author : Fabien Georget
 
 Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
     * Neither the name of the Princeton University nor the
       names of its contributors may be used to endorse or promote products
       derived from this software without specific prior written permission.
 
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 ---------------------------------------------------------*/
 
 #include "diffusion_numbering.hpp"
 #include "reactmicp/meshes/mesh1d.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace diffusion {
 
 
 // Boundary Condition
 
 void BoundaryCondition::fix_all()
 {
     if (not composition.is_valid())
         throw std::runtime_error("The composition has not been correctly initialized");
 
     std::shared_ptr<database::DataContainer> data = composition.get_database();
     component.reserve(data->nb_component - 1);
     tot_mobile_concentration.reserve(data->nb_component -1);
     mineral.reserve(data->nb_mineral);
     for (int comp: data->range_aqueous_component())
     {
         component.push_back(comp);
         tot_mobile_concentration.push_back(comp);
     }
     for (int min: data->range_mineral())
     {
         mineral.push_back(min);
     }
 }
 
 // Diffusion numbering
 
 DiffusionNumbering::DiffusionNumbering(ind_t number_nodes,
                    ind_t number_aqueous_component,
                    ind_t number_mineral, const ListBoundaryCondition &list_bc, Variables &variable):
     nb_aq_component(number_aqueous_component), nb_mineral(number_mineral),
     nb_ndf(2*number_aqueous_component+number_mineral),
     m_ideq(nb_ndf, number_nodes)
 {
     number_equations(list_bc, variable);
 }
 
 
 void DiffusionNumbering::number_equations(const ListBoundaryCondition& list_bc,
                       Variables& variable)
 {
     initialize_variable(variable);
     // First we parse the bcs, to find the correct bc to apply
     parse_bcs(list_bc, variable);
     // Then we number the equations
     do_numbering();
 }
 
 
 void DiffusionNumbering::initialize_variable(Variables& variable)
 {
     variable.displacement = Vector::Zero(m_ideq.rows()*m_ideq.cols());
     variable.velocity = Vector::Zero(m_ideq.rows()*m_ideq.cols());
 }
 
 void DiffusionNumbering::parse_bcs(const ListBoundaryCondition& list_bc,
                                   Variables& variable)
 {
+    variable.displacement.setZero();
+    variable.velocity.setZero();
     m_ideq.setZero();
     m_nodes_with_bc.reserve(list_bc.size());
     for (BoundaryCondition bc: list_bc)
     {
         const ind_t node = bc.node;
         m_nodes_with_bc.push_back(node);
         for (ind_t component: bc.tot_mobile_concentration)
         {
             m_ideq(get_ndof_diffusion(component), node) = no_equation;
             mobile_total_concentration(node, component, variable) = bc.composition.total_aqueous_concentration_component(component);
         }
         for (ind_t component: bc.component)
         {
             m_ideq(get_ndof_massbalance(component), node) = no_equation;
             component_concentration(node, component, variable) = std::log10(bc.composition.molality_component(component));
         }
         for (ind_t mineral: bc.mineral)
         {
             m_ideq(get_ndof_mineral(mineral), node) = no_equation;
             mole_mineral(node, mineral, variable) = bc.composition.moles_mineral(mineral);
         }
     }
+
+    //
+    for (int node=0; node<m_ideq.cols(); ++node)
+    {
+        for (ind_t component=1; component<5; ++component)
+        {
+        //    m_ideq(get_ndof_diffusion(component), node) = no_equation;
+//            m_ideq(get_ndof_massbalance(component), node) = no_equation;
+        }
+//        for (ind_t mineral=0; mineral<5; ++mineral)
+//        {
+//            m_ideq(get_ndof_mineral(mineral), node) = no_equation;
+//        }
+    }
 }
 
 void DiffusionNumbering::do_numbering()
 {
     ind_t neq = 0;
     for (ind_t node=0; node<m_ideq.cols(); ++node)
     {
         for (ind_t df=0; df<get_ndf(); ++df)
         {
             if (m_ideq(df, node) != no_equation)
             {
                 m_ideq(df, node) = neq;
                 ++neq;
             }
         }
     }
     nb_neq = neq;
+    std::cout << m_ideq << std::endl;
 }
 
 
 //! \brief Return true if the node has a boundary condition
 double DiffusionNumbering::node_has_bc(ind_t node)
 {
     auto search = std::find(m_nodes_with_bc.begin(), m_nodes_with_bc.end(), node);
     if (search == m_nodes_with_bc.end()) return false;
     else return true;
 }
 
 } // end namespace diffusion
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
diff --git a/src/reactmicp/systems/diffusion/diffusion_numbering.hpp b/src/reactmicp/systems/diffusion/diffusion_numbering.hpp
index 2a28e60..e7232e8 100644
--- a/src/reactmicp/systems/diffusion/diffusion_numbering.hpp
+++ b/src/reactmicp/systems/diffusion/diffusion_numbering.hpp
@@ -1,202 +1,202 @@
 /*-------------------------------------------------------
 
  - Module : reactmicp/systems/diffusion
  - File : diffusion_numbering.hpp
  - Author : Fabien Georget
 
 Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
     * Neither the name of the Princeton University nor the
       names of its contributors may be used to endorse or promote products
       derived from this software without specific prior written permission.
 
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 ---------------------------------------------------------*/
 
 #ifndef SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_NUMBERING_HPP
 #define SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_NUMBERING_HPP
 
 
 #include <memory>
 #include "reactmicp/common.hpp"
 #include "specmicp/equilibrium_data.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 
 namespace mesh {
 // forward declaration
 class Mesh1D;
 } // end namespace mesh
 
 namespace systems {
 
 namespace diffusion {
 
 struct BoundaryCondition
 {
     int node;
     EquilibriumState composition;
     std::vector<int> tot_mobile_concentration;
     std::vector<int> component;
     std::vector<int> mineral;
 
     void fix_all();
 };
 
 using ListBoundaryCondition = std::vector<BoundaryCondition>;
 
 //! \brief Equation numbering scheme
 class DiffusionNumbering
 {
 public:
     const int no_equation = -1;
     using Variables = ParabolicVariables;
 
     DiffusionNumbering(ind_t number_nodes,
                        std::shared_ptr<database::DataContainer> data,
                        const ListBoundaryCondition& list_bc,
                        Variables& variable):
         DiffusionNumbering(number_nodes, data->nb_component-1, data->nb_mineral, list_bc, variable)
     {}
 
     DiffusionNumbering(ind_t number_nodes,
                        ind_t number_aqueous_component,
                        ind_t number_mineral,
                        const ListBoundaryCondition& list_bc,
                        Variables& variable);
 
     //! \brief Return the number of equations
     ind_t get_neq() const {return nb_neq;}
     //! \brief Return the number of degree of freedom (per node)
     ind_t get_ndf() const {return nb_ndf;}
 
     // Degree of freedom numbering
     // ===========================
 
     //! \brief Return the number of the dof
     ind_t get_dof(ind_t node,  ind_t elemental_dof) const {
         return elemental_dof + nb_ndf*node;}
 
     //! \brief Return the ndof of a diffusion equation
     ind_t get_ndof_diffusion(ind_t ind_component) const {
         return ind_component-1;}
 
     //! \brief Return the ndof of a mass balance equation
     ind_t get_ndof_massbalance(ind_t ind_component) const {
         return nb_aq_component + ind_component - 1;}
 
     //! \brief Return the ndof of a mineral equation
     ind_t get_ndof_mineral(ind_t ind_mineral) const {
         return 2*nb_aq_component + ind_mineral;}
 
     //! \brief Return the dof of a diffusion equation
     ind_t get_dof_diffusion(ind_t node, ind_t ind_component) const {
         return get_dof(node, get_ndof_diffusion(ind_component));}
 
     //! \brief Return the dof of a mass balance equation
     ind_t get_dof_massbalance(ind_t node, ind_t ind_component) const {
         return get_dof(node, get_ndof_massbalance(ind_component));}
 
     //! \brief Return the dof of a mineral equation
     ind_t get_dof_mineral(ind_t node, ind_t ind_mineral) const {
         return get_dof(node, get_ndof_mineral(ind_mineral));}
 
     //! \brief Return the index of the equation
     //!
     //! Return no_equation if there is no equation for this degree of freedom
     ind_t id_equation(ind_t node, ind_t elemental_dof) const {return m_ideq(elemental_dof, node);}
 
     //! \brief Return the equation's index of a diffusion equation
     ind_t id_equation_diffusion(ind_t node, ind_t ind_component) const {
         return id_equation(node, get_ndof_diffusion(ind_component));}
 
     //! \brief Return the equation's index of a mass balance equation
     ind_t id_equation_massbalance(ind_t node, ind_t ind_component) const {
         return id_equation(node, get_ndof_massbalance(ind_component));}
 
     //! \brief Return the equation's index of a mineral equation
     ind_t id_equation_mineral(ind_t node, ind_t ind_mineral) const {
         return id_equation(node, get_ndof_mineral(ind_mineral));}
 
     // Access primary variable
     // =======================
 
     //! \brief Return the total mobile concentration for 'component' at 'node'
     double mobile_total_concentration(ind_t node, ind_t component, const Variables& var) const {
         return var.displacement(get_dof_diffusion(node, component));
     }
     //! \brief Return a reference to the total mobile concentration for 'component' at 'node'
     double& mobile_total_concentration(ind_t node, ind_t component, Variables& var) {
         return var.displacement.coeffRef(get_dof_diffusion(node, component));
     }
 
 
     //! \brief Return the component concentration for 'component' at 'node'
     double component_concentration(ind_t node, ind_t component, const Variables& var) const {
         return var.displacement(get_dof_massbalance(node, component));
     }
     //! \brief Return a reference to the component concentration for 'component' at 'node'
     double& component_concentration(ind_t node, ind_t component, Variables& var) {
-        return var.displacement(get_dof_massbalance(node, component));
+        return var.displacement.coeffRef(get_dof_massbalance(node, component));
     }
 
     //! \brief Return the mineral mole number for 'mineral' at 'node'
     double mole_mineral(ind_t node, ind_t mineral, const Variables& var) const {
         return var.displacement(get_dof_mineral(node, mineral));
     }
     //! \brief Return a reference to the mineral mole number for 'mineral' at 'node'
     double& mole_mineral(ind_t node, ind_t mineral, Variables& var) {
         return var.displacement.coeffRef(get_dof_mineral(node, mineral));
     }
 
     //! \brief Return true if the node has a boundary condition
     double node_has_bc(ind_t node);
 
 private:
     //! \brief Initialize the variable
     void initialize_variable(Variables& variable);
 
     //! \brief Parse the BCs and number the equations
     void number_equations(const ListBoundaryCondition& list_bc,
                           Variables& variable);
     //! \brief Parse the BCS, set the corresponding values in variable
     void parse_bcs(const ListBoundaryCondition& list_bc,
                    Variables& variable);
     //! \brief number the equations
     void do_numbering();
 
     ind_t nb_aq_component;
     ind_t nb_mineral;
 
     ind_t nb_ndf;
     ind_t nb_neq;
 
     Eigen::MatrixXi m_ideq;
     std::vector<ind_t> m_nodes_with_bc;
 };
 
 } // end namespace diffusion
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
 
 #endif // SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_NUMBERING_HPP
diff --git a/src/reactmicp/systems/diffusion/diffusion_secondary.cpp b/src/reactmicp/systems/diffusion/diffusion_secondary.cpp
index 5a77d66..430bffc 100644
--- a/src/reactmicp/systems/diffusion/diffusion_secondary.cpp
+++ b/src/reactmicp/systems/diffusion/diffusion_secondary.cpp
@@ -1,173 +1,181 @@
 /*-------------------------------------------------------
 
  - Module : reactmicp/systems/diffusion
  - File : diffusion_secondary.cpp
  - Author : Fabien Georget
 
 Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
     * Neither the name of the Princeton University nor the
       names of its contributors may be used to endorse or promote products
       derived from this software without specific prior written permission.
 
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 ---------------------------------------------------------*/
 
 #include "diffusion_secondary.hpp"
 #include "physics/laws.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace diffusion {
 
 
 // ================================== //
 //                                    //
 //  Secondary variables wibbly-timey  //
 //                                    //
 // ================================== //
 
 //! \brief Set secondary concentrations
 void DiffusionSecondaryVariables::compute_secondary_concentrations(const Variables& variable)
 {
     for (ind_t node=0; node<m_secondary_concentration.cols(); ++node)
     {
         nodal_secondary_concentrations(node, variable);
     }
 }
 
 //! \brief Set secondary concentrations for node 'node'
 void DiffusionSecondaryVariables::nodal_secondary_concentrations(ind_t node, const Variables& variable)
 {
     for (int species: m_data->range_aqueous())
     {
         double sum = -m_data->logk_aqueous(species) - loggamma_secondary(node, species);
         for (int component: m_data->range_aqueous_component())
         {
             sum +=  m_data->nu_aqueous(species, component)*(
-                        loggamma_component(node,component)
+                        loggamma_component(node, component)
                         + m_ideq.component_concentration(node, component, variable)
                     );
         }
         secondary_concentration(node, species) = pow10(sum);
     }
 }
 
+double DiffusionSecondaryVariables::norm_secondary_variables()
+{
+    return  m_secondary_variables.block(
+                1, 0,
+                m_data->nb_aqueous+m_data->nb_component,
+                m_secondary_concentration.cols()
+            ).norm();
+}
+
+
 //! \brief Solve for secondary variables
 int DiffusionSecondaryVariables::solve_secondary_variables(const Variables& variable)
 {
     bool may_have_converged = false;
-    double previous_norm = m_secondary_variables.block(
-                0, 0, m_data->nb_aqueous+m_data->nb_component, m_secondary_concentration.cols()).norm();
+    double previous_norm = norm_secondary_variables();
     for (int i=0; i<6; ++i)
     {
         compute_secondary_concentrations(variable);
         compute_secondary_variables(variable);
-        double new_norm = m_secondary_variables.block(
-                    0, 0, m_data->nb_aqueous+m_data->nb_component, m_secondary_concentration.cols()).norm();
-        if (std::abs(previous_norm - new_norm) < 1e-6) {may_have_converged=true; break;}
+        double new_norm = norm_secondary_variables();
+        if (std::abs(previous_norm - new_norm)/previous_norm < 1e-6) {may_have_converged=true; break;}
         previous_norm = new_norm;
     }
     return may_have_converged;
 }
 
 //! \brief Compute secondary variables
 void DiffusionSecondaryVariables::compute_secondary_variables(const Variables& variable)
 {
     compute_ionic_strength(variable);
     compute_loggamma();
 }
 
 void DiffusionSecondaryVariables::compute_ionic_strength(const Variables& variable)
 {
     for (ind_t node=0; node<m_secondary_concentration.cols(); ++node)
     {
         nodal_ionic_strength(node, variable);
     }
 }
 
 //! \brief compute the ionic strength for node 'node'
 void DiffusionSecondaryVariables::nodal_ionic_strength(ind_t node, const Variables& variable)
 {
     double ionics = 0;
     for (int component: m_data->range_aqueous_component())
     {
         if (m_data->param_aq(component, 0) == 0) continue;
         ionics += std::pow(m_data->param_aq(component, 0), 2)*pow10(
                     m_ideq.component_concentration(node, component, variable));
     }
     for (int aqueous: m_data->range_aqueous())
     {
         ind_t idaq = m_data->nb_component + aqueous;
         if (m_data->param_aq(idaq, 0) == 0) continue;
         ionics += std::pow(m_data->param_aq(idaq, 0), 2)*secondary_concentration(node, aqueous);
     }
      ionic_strength(node) = ionics/2;
 }
 
 void DiffusionSecondaryVariables::compute_loggamma()
 {
     for (ind_t node=0; node<m_secondary_concentration.cols(); ++node)
     {
         nodal_loggamma(node);
     }
 }
 
 //! \brief compute the logarithm of the activity coefficients for node 'node'
 void DiffusionSecondaryVariables::nodal_loggamma(ind_t node)
 {
     nodal_loggamma_component(node);
     nodal_loggamma_aqueous(node);
 }
 
 
 //! \brief compute the logarithm of the activity coefficients for node 'node'
 void DiffusionSecondaryVariables::nodal_loggamma_component(ind_t node)
 {
-    double sqrti = std::sqrt(ionic_strength(node));
+    const double sqrti = std::sqrt(ionic_strength(node));
     for (int component: m_data->range_aqueous_component())
     {
         loggamma_component(node, component) = laws::extended_debye_huckel(ionic_strength(node), sqrti,
                                                                       m_data->param_aq(component, 0),
                                                                       m_data->param_aq(component, 1),
                                                                       m_data->param_aq(component, 2));
     }
 }
 
 //! \brief compute the logarithm of the activity coefficients for node 'node'
 void DiffusionSecondaryVariables::nodal_loggamma_aqueous(ind_t node)
 {
-    double sqrti = std::sqrt(ionic_strength(node));
+    const double sqrti = std::sqrt(ionic_strength(node));
     for (int aqueous: m_data->range_aqueous())
     {
         ind_t idaq = m_data->nb_component+aqueous;
         loggamma_secondary(node, aqueous) = laws::extended_debye_huckel(ionic_strength(node), sqrti,
                                                                         m_data->param_aq(idaq, 0),
                                                                         m_data->param_aq(idaq, 1),
                                                                         m_data->param_aq(idaq, 2));
     }
 }
 
 } // end namespace diffusion
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
diff --git a/src/reactmicp/systems/diffusion/diffusion_secondary.hpp b/src/reactmicp/systems/diffusion/diffusion_secondary.hpp
index a2e7e3f..088db65 100644
--- a/src/reactmicp/systems/diffusion/diffusion_secondary.hpp
+++ b/src/reactmicp/systems/diffusion/diffusion_secondary.hpp
@@ -1,161 +1,173 @@
 /*-------------------------------------------------------
 
  - Module : reactmicp/systems/diffusion
  - File : diffusion_secondary.hpp
  - Author : Fabien Georget
 
 Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
     * Neither the name of the Princeton University nor the
       names of its contributors may be used to endorse or promote products
       derived from this software without specific prior written permission.
 
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 ---------------------------------------------------------*/
 
 #ifndef SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_SECONDARY_HPP
 #define SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_SECONDARY_HPP
 
 
 #include <memory>
 
 #include "database/data_container.hpp"
 
 #include "reactmicp/common.hpp"
 #include "diffusion_numbering.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace diffusion {
 
 //! \brief handle secondary variables for the diffusion system
 class DiffusionSecondaryVariables
 {
 public:
     using Variables = ParabolicVariables;
 
     DiffusionSecondaryVariables(int nb_nodes,
                                 std::shared_ptr<database::DataContainer> thedatabase,
                                 DiffusionNumbering& numberer):
         m_data(thedatabase), m_ideq(numberer),
         m_secondary_variables(Eigen::MatrixXd(thedatabase->nb_component+thedatabase->nb_aqueous+1, nb_nodes)),
         m_secondary_concentration(Eigen::MatrixXd(thedatabase->nb_aqueous, nb_nodes))
     {
         m_secondary_variables.setZero();
         m_secondary_concentration.setZero();
     }
 
     // Secondary variables
     // ===================
 
     // Algorithms
     // ----------
 
     //! \brief Set secondary concentrations
     void compute_secondary_concentrations(const Variables& variable);
 
     //! \brief Set secondary concentrations for node 'node'
     void nodal_secondary_concentrations(ind_t node, const Variables& variable);
 
     //! \brief Solve for secondary variables
     int solve_secondary_variables(const Variables& variable);
 
     //! \brief Compute secondary variables
     void compute_secondary_variables(const Variables& variable);
 
     //! \brief compute the ionic strengths
     void compute_ionic_strength(const Variables& variable);
 
     //! \brief compute the ionic strength for node 'node'
     void nodal_ionic_strength(ind_t node, const Variables& variable);
 
     //! \brief compute the activity coefficients
     void compute_loggamma();
 
     //! \brief compute the logarithm of the activity coefficients for node 'node'
     void nodal_loggamma(ind_t node);
 
     //! \brief compute the activity coefficient for the components of node 'node'
     void nodal_loggamma_component(ind_t node);
 
     //! \brief compute the activity coefficient for the aqueous species of node 'node'
     void nodal_loggamma_aqueous(ind_t node);
 
     // Getters and setters
     // -------------------
 
     //! \brief Secondary species concentration
     double& secondary_concentration(ind_t node, ind_t aqueous) {
+        assert(aqueous >= 0 and aqueous < m_data->nb_aqueous);
         return m_secondary_concentration.coeffRef(aqueous, node);
     }
     //! \brief Secondary species concentration
     double secondary_concentration(ind_t node, ind_t aqueous) const {
+        assert(aqueous >= 0 and aqueous < m_data->nb_aqueous);
         return m_secondary_concentration.coeff(aqueous, node);
     }
 
     //! \brief Return the loggamma of a component
     double& loggamma_component(ind_t node, ind_t component) {
+        assert(component >= 1 and component < m_data->nb_component);
         return m_secondary_variables.coeffRef(component, node);
     }
-
     //! \brief Return the loggamma of a component
     double loggamma_component(ind_t node, ind_t component) const {
+        assert(component >= 1 and component < m_data->nb_component);
         return m_secondary_variables.coeff(component, node);
     }
 
     //! \brief Return the loggamma of a secondary species
-    double& loggamma_secondary(ind_t node, ind_t species) {
-        return m_secondary_variables.coeffRef(m_data->nb_component + species, node);
+    double& loggamma_secondary(ind_t node, ind_t aqueous) {
+        assert(aqueous >= 0 and aqueous < m_data->nb_aqueous);
+        return m_secondary_variables.coeffRef(m_data->nb_component + aqueous, node);
     }
-
     //! \brief Return the loggamma of a secondary species
-    double loggamma_secondary(ind_t node, ind_t species) const {
-        return m_secondary_variables.coeff(m_data->nb_component + species, node);
+    double loggamma_secondary(ind_t node, ind_t aqueous) const {
+        assert(aqueous >= 0 and aqueous < m_data->nb_aqueous);
+        return m_secondary_variables.coeff(m_data->nb_component + aqueous, node);
     }
 
     //! \brief Return ionic strength
     double& ionic_strength(ind_t node) {
         return m_secondary_variables.coeffRef(m_data->nb_component+m_data->nb_aqueous, node);
     }
     //! \brief Return ionic strength
     double ionic_strength(ind_t node) const {
         return m_secondary_variables.coeff(m_data->nb_component+m_data->nb_aqueous, node);
     }
 
+
+    //! \brief Return norm of the loggamma
+    double loggamma_norm() const {
+        return m_secondary_variables.norm();
+    }
+
 private:
+    // Return the norm used to check the convergence of the fixed-point iterations
+    double norm_secondary_variables();
 
     std::shared_ptr<database::DataContainer> m_data;
     DiffusionNumbering& m_ideq;
 
     Eigen::MatrixXd m_secondary_variables;
     Eigen::MatrixXd m_secondary_concentration;
 
 };
 
 
 } // end namespace diffusion
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
 
 #endif // SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_SECONDARY_HPP
 
diff --git a/tests/reactmip/systems/diffusion/solving.cpp b/tests/reactmip/systems/diffusion/solving.cpp
index a101197..45e0b6f 100644
--- a/tests/reactmip/systems/diffusion/solving.cpp
+++ b/tests/reactmip/systems/diffusion/solving.cpp
@@ -1,44 +1,68 @@
 #include "catch.hpp"
 
 #include "utils.hpp"
 #include "reactmicp/systems/diffusion/diffusion.hpp"
 #include "reactmicp/meshes/mesh1d.hpp"
 
 #include "database/database.hpp"
 
 #include "reactmicp/micpsolvers/micpsolver.hpp"
 #include "utils/log.hpp"
 #include <iostream>
 using namespace specmicp;
 using namespace specmicp::reactmicp;
 using namespace specmicp::reactmicp::systems;
 using namespace specmicp::reactmicp::systems::diffusion;
 
 TEST_CASE("Solving the diffusion system", "[Diffusion,MiCP]") {
     specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
     specmicp::logger::ErrFile::stream() = &std::cerr;
 
     DiffusionSecondaryVariables::Variables var;
     std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
     EquilibriumState comp_bc;
     comp_bc = blank_composition(database);
     BoundaryCondition bc_node0;
     bc_node0.node = 0;
     bc_node0.composition = comp_bc;
     bc_node0.fix_all();
     ListBoundaryCondition bcs = {bc_node0,};
 
     EquilibriumState initial_state = sample_carbo_composition(database);
 
     std::shared_ptr<mesh::Mesh1D> themesh = std::make_shared<mesh::Mesh1D>(4, 1e-3, 0.001);
     std::shared_ptr<DiffusionParameter> parameters = std::make_shared<DiffusionParameter>(1e-8, 0.2);
 
     SECTION("Try solving") {
-        auto program = std::make_shared<DiffusionProgram>(themesh, database, parameters, bcs, var);
+        auto program = std::make_shared<DiffusionProgram>(themesh, database, parameters, bcs, initial_state, var);
 
         specmicp::reactmicp::micpsolver::MiCPParabolicSolver<DiffusionProgram> solver(program);
 
-        auto retcode = solver.solve_one_timestep(0.0001, var);
-        CHECK((int) retcode == 2);
+
+        int ndf = program->get_ndf();
+        for (int i=0; i<program->get_ndf(); ++i)
+        {
+        for (int node: themesh->range_nodes())
+        {
+            std::cout << var.displacement(i+ndf*node) << "\t";
+        }
+        std::cout << std::endl;
+        }
+        //std::cout << var.displacement << std::endl;
+
+        for (int k=0; k<10; ++k)
+        {
+        auto retcode = solver.solve_one_timestep(2000, var);
+
+        for (int i=0; i<program->get_ndf(); ++i)
+        {
+        for (int node: themesh->range_nodes())
+        {
+            std::cout << var.displacement(i+ndf*node) << "\t";
+        }
+        std::cout << std::endl;
+        }
+        REQUIRE((int) retcode > 0);
+        }
     }
 }