diff --git a/src/dfpm/1dtransport/diffusion.cpp b/src/dfpm/1dtransport/diffusion.cpp
index 93b2c66..aedac2f 100644
--- a/src/dfpm/1dtransport/diffusion.cpp
+++ b/src/dfpm/1dtransport/diffusion.cpp
@@ -1,179 +1,194 @@
 #include "diffusion.hpp"
 #include "diffusion_parameters.hpp"
 #include "dfpm/meshes/mesh1d.hpp"
-#include "dfpm/types.hpp"
 
 namespace specmicp {
 namespace dfpm {
 
 // SaturatedDiffusion1D::
 
 SaturatedDiffusion1D::SaturatedDiffusion1D(
         mesh::Mesh1DPtr the_mesh,
-        SaturatedDiffusion1DParameters parameters,
+        std::shared_ptr<SaturatedDiffusion1DParameters> parameters,
         std::vector<index_t> list_bcs
         ):
-    m_ndf(1),
     m_tot_ndf(the_mesh->nb_nodes()),
     m_mesh(the_mesh),
     m_param(parameters),
+    m_internal_flow(Vector::Zero(m_tot_ndf)),
     m_external_flow(Vector::Zero(m_tot_ndf))
 {
     number_equations(list_bcs);
 }
 
 void SaturatedDiffusion1D::number_equations(std::vector<index_t> list_bcs)
 {
-    m_id_equations = Eigen::VectorXi::Zeros(m_tot_ndf);
-    for (std::vector<index_t>::Iterator it=list_bcs.begin(); it!=list_bcs.end(); ++it)
+    m_id_equations = Eigen::VectorXi::Zero(m_tot_ndf);
+    for (auto it=list_bcs.begin(); it!=list_bcs.end(); ++it)
     {
         m_id_equations(*it) = no_equation;
     }
     index_t neq = 0;
     for (index_t node: m_mesh->range_nodes())
     {
         if (m_id_equations(node) == no_equation) continue;
         m_id_equations(node) = neq;
         ++ neq;
     }
     m_neq = neq;
 }
 
 void SaturatedDiffusion1D::compute_residuals(
         const Vector& displacement,
         const Vector& velocity,
         Vector& residual
         )
 {
+    m_internal_flow.setZero();
     residual.resize(get_neq());
     residual.setZero();
     for (index_t element: m_mesh->range_elements())
     {
-        Vector element_residuals(Vector::Eigen::Zero(2));
-        element_residuals(displacement, velocity, element_residuals);
+        Vector elem_residuals(2);
+        elem_residuals.setZero();
+        element_residuals(element, displacement, velocity, elem_residuals);
         for (index_t enode: m_mesh->range_nen())
         {
             if (id_equation(m_mesh->get_node(element, enode)) == no_equation) continue;
-            residual(id_equation(m_mesh->get_node(element, enode))) += element_residuals(enode);
+            residual(id_equation(m_mesh->get_node(element, enode))) += elem_residuals(enode);
         }
     }
 }
 
-
-
 void SaturatedDiffusion1D::element_residuals(
+        index_t element,
         const Vector& displacement,
         const Vector& velocity,
-        Vector& element_residual
+        Vector& elem_residuals
         )
 {
     Eigen::Matrix<scalar_t, 2, 2> jacob;
     Eigen::Matrix<scalar_t, 2, 1> evelocity, edisplacement;
     scalar_t mass_coeff = -(m_mesh->get_volume_element(element)/2.0);
 
     const index_t node_0 = m_mesh->get_node(element, 0);
     const index_t node_1 = m_mesh->get_node(element, 1);
 
     const scalar_t porosity = (m_param->porosity(node_0) +
                                m_param->porosity(node_1) )/2.0;
     const scalar_t diff_coeff = 1.0/(0.5/m_param->diffusion_coefficient(node_0) +
                                      0.5/m_param->diffusion_coefficient(node_1));
 
     const index_t dof_0 = node_0;
     const index_t dof_1 = node_1;
 
     scalar_t flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element)
-                             * porosity
+                            // * porosity
                              * diff_coeff
                              );
     jacob << 1.0, -1.0, -1.0, 1.0;
     jacob *= flux_coeff;
     evelocity << velocity(dof_0)*
                  mass_coeff*m_param->porosity(node_0),
             velocity(dof_1)*
             mass_coeff*m_param->porosity(node_1);
 
     edisplacement << displacement(dof_0),
             displacement(dof_1);
 
-    element_residual = jacob*edisplacement;
+    elem_residuals = jacob*edisplacement;
 
-    m_internal_flow(dof_0) += element_residual(0);
-    m_internal_flow(dof_1) += element_residual(1);
+    m_internal_flow(dof_0) += elem_residuals(0);
+    m_internal_flow(dof_1) += elem_residuals(1);
 
     for (index_t en: m_mesh->range_nen())
     {
-        element_residual(en) += evelocity(en);
-        element_residual(en) += (m_mesh->get_volume_element(element)/2.0
-                                 *external_flow(m_mesh->get_node(element, en), component));
+        elem_residuals(en) += evelocity(en);
+        elem_residuals(en) += (m_mesh->get_volume_element(element)/2.0
+                                 *external_flow(m_mesh->get_node(element, en)));
     }
 }
 
 void SaturatedDiffusion1D::compute_jacobian(
         Vector& displacement,
         Vector& velocity,
         Eigen::SparseMatrix<scalar_t>& jacobian,
         scalar_t alphadt
         )
 {
 
     list_triplet_t jacob;
     const index_t estimation = m_mesh->nb_nodes()*(m_mesh->nen);
     jacob.reserve(estimation);
     for (index_t element: m_mesh->range_elements())
     {
-        element_jacobian_transport(element, displacement, velocity, jacob, alphadt);
+        element_jacobian(element, displacement, velocity, jacob, alphadt);
     }
     jacobian = Eigen::SparseMatrix<scalar_t>(get_neq(), get_neq());
     jacobian.setFromTriplets(jacob.begin(), jacob.end());
 }
 
-void SaturatedDiffusion1D::element_jacobian_transport(
+void SaturatedDiffusion1D::element_jacobian(
         index_t element,
         Vector& displacement,
         Vector& velocity,
         list_triplet_t& jacobian,
         scalar_t alphadt)
 {
     Eigen::VectorXd element_residual_orig(Eigen::VectorXd::Zero(2));
     element_residuals(element, displacement, velocity, element_residual_orig);
 
     for (index_t en: m_mesh->range_nen())
     {
         Eigen::VectorXd element_residual(Eigen::VectorXd::Zero(2));
 
         const index_t node = m_mesh->get_node(element, en);
         const index_t dof = node;
         const index_t idc = id_equation(dof);
 
         if (idc == no_equation) continue;
 
         const scalar_t tmp_v = velocity(dof);
         const scalar_t tmp_d = displacement(dof);
 
         scalar_t h = eps_jacobian*std::abs(tmp_v);
         if (h < 1e-4*eps_jacobian) h = eps_jacobian;
         velocity(dof) = tmp_v + h;
         h = velocity(dof) - tmp_v;
 
         displacement(dof) = tmp_d + alphadt*h;
 
         element_residuals(element, displacement, velocity, element_residual);
         velocity(dof) = tmp_v;
         displacement(dof) = tmp_d;
 
         for (index_t enr: m_mesh->range_nen())
         {
             const index_t noder = m_mesh->get_node(element, enr);
             const index_t idr = id_equation(noder);
 
             if (idr == no_equation) continue;
             jacobian.push_back(triplet_t(idr, idc, (element_residual(enr) - element_residual_orig(enr))/h));
         }
     }
+}
 
-
+void SaturatedDiffusion1D::update_solution(const Vector &update,
+        scalar_t lambda,
+        scalar_t alpha_dt,
+        Vector &predictor,
+        Vector &displacement,
+        Vector &velocity
+        )
+{
+    for (index_t node: m_mesh->range_nodes())
+    {
+        const index_t id = id_equation(node);
+        if (id == no_equation) continue;
+        velocity(node) += lambda*update(id);
+    }
+    displacement = predictor + alpha_dt*velocity;
 }
 
 } // end namespace dfpm
 } // end namespace specmicp
diff --git a/src/dfpm/1dtransport/diffusion.hpp b/src/dfpm/1dtransport/diffusion.hpp
index 2532de7..6871df4 100644
--- a/src/dfpm/1dtransport/diffusion.hpp
+++ b/src/dfpm/1dtransport/diffusion.hpp
@@ -1,92 +1,108 @@
 #ifndef SPECMICP_DFPM_1DTRANSPORT_DIFFUSION_HPP
 #define SPECMICP_DFPM_1DTRANSPORT_DIFFUSION_HPP
 
 #include <memory>
 
 #include "common.hpp"
+
+#include "dfpm/types.hpp"
 #include "dfpmsolver/parabolic_program.hpp"
 #include "dfpm/meshes/mesh1dfwd.hpp"
 
 namespace specmicp {
 
 namespace dfpm {
 
 class SaturatedDiffusion1DParameters;
 
 class SaturatedDiffusion1D: public dfpmsolver::ParabolicProgram<SaturatedDiffusion1D>
 {
 public:
+    SaturatedDiffusion1D(
+            mesh::Mesh1DPtr the_mesh,
+            std::shared_ptr<SaturatedDiffusion1DParameters>  parameters,
+            std::vector<index_t> list_bcs
+            );
+
     //! \brief Return the number of equations
     index_t get_neq() const {return m_neq;}
 
     //! \brief Return the number of degrees of freedom per node
     index_t get_ndf() const {return 1;}
 
     //! \brief Return the total number of degrees of freedom
     index_t get_tot_ndf() const {return m_tot_ndf;}
 
     //! \brief Return the id of the equation corresponding to the degree of freedom 'id_dof'
     //!
     //! Return 'no_equation' if no equation exist
     index_t id_equation(index_t id_dof) const {return m_id_equations(id_dof);}
 
-    void element_residuals(const Vector& displacement,
+    void element_residuals(index_t element,
+                           const Vector& displacement,
                            const Vector& velocity,
                            Vector& element_residual
                            );
 
     //! \brief Compute the residuals
     void compute_residuals(const Vector& displacement,
                            const Vector& velocity,
                            Vector& residual
                            );
 
     //! \brief Compute the jacobian
     void compute_jacobian(Vector& displacement,
                           Vector& velocity,
                           Eigen::SparseMatrix<scalar_t>& jacobian,
                           scalar_t alphadt
                           );
+    void element_jacobian(
+            index_t element,
+            Vector& displacement,
+            Vector& velocity,
+            list_triplet_t& jacobian,
+            scalar_t alphadt);
 
     //! \brief Update the solutions
     void update_solution(const Vector& update,
                          scalar_t lambda,
                          scalar_t alpha_dt,
                          Vector& predictor,
                          Vector& displacement,
                          Vector& velocity);
 
     //! \brief Apply boundary conditions to the velocity vector
     //!
     //! by default do nothing.
     void apply_bc(scalar_t dt,
                   const Vector& displacement,
                   Vector& velocity)
     {}
 
     //! \brief Return the value of the external flow for dof 'id_dof'
     scalar_t external_flow(index_t id_dof) const {return m_external_flow(id_dof);}
     //! \brief Return a reference to the value of the external flow for dof 'id_dof'
     scalar_t& external_flow(index_t id_dof) {return m_external_flow(id_dof);}
     //! \brief Return a reference to the vector of external flow
     Vector& external_flow() {return m_external_flow;}
 
 private:
     void number_equations(std::vector<index_t> list_bcs);
 
     index_t m_tot_ndf;
     index_t m_neq;
     Eigen::VectorXi m_id_equations;
 
     mesh::Mesh1DPtr m_mesh;
     std::shared_ptr<SaturatedDiffusion1DParameters> m_param;
 
+    Vector m_internal_flow;
     Vector m_external_flow;
 };
 
 
 } // end namespace dfpm
 
 } // end namespace specmicp
 
 #endif // SPECMICP_DFPM_1DTRANSPORT_DIFFUSION_HPP
diff --git a/src/dfpm/1dtransport/diffusion_parameters.hpp b/src/dfpm/1dtransport/diffusion_parameters.hpp
index 8478551..5d81787 100644
--- a/src/dfpm/1dtransport/diffusion_parameters.hpp
+++ b/src/dfpm/1dtransport/diffusion_parameters.hpp
@@ -1,26 +1,25 @@
 #ifndef SPECMICP_DFPM_1DTRANSPORT_DIFFUSIONPARAMETERS_HPP
 #define SPECMICP_DFPM_1DTRANSPORT_DIFFUSIONPARAMETERS_HPP
 
 #include "common.hpp"
 
 namespace specmicp {
 namespace dfpm {
 
 struct SaturatedDiffusion1DParameters
 {
     Vector porosity;
-    Vector diffusion_coefficients;
-    scalar_t density_water;
+    Vector diffusion_coefficient;
 
-    SaturatedDiffusion1DParameters(nb_nodes):
+    SaturatedDiffusion1DParameters(index_t nb_nodes):
         porosity(nb_nodes),
-        diffusion_coefficients(nb_nodes)
+        diffusion_coefficient(nb_nodes)
     {}
 
 
 };
 
 } // end namespace dfpm
 } // end namespace specmicp
 
 #endif // SPECMICP_DFPM_1DTRANSPORT_DIFFUSIONPARAMETERS_HPP