diff --git a/python/tamaas/nonlinear_solvers/scipy_solver.py b/python/tamaas/nonlinear_solvers/scipy_solver.py
index 9cdcbe0..6b6930c 100644
--- a/python/tamaas/nonlinear_solvers/scipy_solver.py
+++ b/python/tamaas/nonlinear_solvers/scipy_solver.py
@@ -1,139 +1,146 @@
 """
 @author Lucas Frérot <lucas.frerot@epfl.ch>
 
 @section LICENSE
 
 Copyright (©)  2018 EPFL  (Ecole Polytechnique  Fédérale de
 Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
 Solides)
 
 Tamaas is free  software: you can redistribute it and/or  modify it under the
 terms  of the  GNU Lesser  General Public  License as  published by  the Free
 Software Foundation, either version 3 of the License, or (at your option) any
 later version.
 
 Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
 details.
 
 You should  have received  a copy  of the GNU  Lesser General  Public License
 along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
 """
 
 from .. import EPSolver
 
 import numpy as np
 from scipy.optimize import newton_krylov, root
 
 __all__ = ["DFSANESolver", "NewtonKrylovSolver", "ToleranceManager"]
 
 
 class ScipySolver(EPSolver):
     """
     Base class for solvers wrapping SciPy routines
     """
 
     def __init__(self, residual, model, callback=None):
         super().__init__(residual)
         self.model = model
         self.callback = callback
         self._x = self.getStrainIncrement()
         self._residual = self.getResidual()
         self.surface_correction = model.getDisplacement()[0, :, :, 2]
         self.last_residual_disp = np.zeros_like(self.surface_correction)
+        self.options = {'ftol': 0, 'fatol': 1e-9}
+
+    @property
+    def tolerance(self):
+        return self.options['fatol']
+
+    @tolerance.setter
+    def set_tolerance(self, v):
+        self.options['fatol'] = v
 
     def solve(self):
         """
         Solve the nonlinear plasticity equation using the scipy_solve routine
         """
         # For initial guess, compute the strain due to boundary tractions
         # self._residual.computeResidual(self._x)
         # self._x[...] = self._residual.getVector()
 
         # Scipy root callback
         def compute_residual(x):
             self._residual.computeResidual(x)
             return self._residual.getVector().copy()
 
         # Solve
         self._x[...] = self.scipy_solve(compute_residual)
 
         # Computing displacements
         self._residual.computeResidualDisplacement(self._x)
 
     def reset(self):
         self.last_residual_disp[...] = 0
         self._x[...] = 0
 
 
 class NewtonKrylovSolver(ScipySolver):
     """
     Solve using a finite-difference Newton-Krylov method
     """
     def __init__(self, residual, model, callback=None):
         ScipySolver.__init__(self, residual, model, callback)
-        self.options = {'ftol': 0, 'fatol': 1e-9}
 
     def scipy_solve(self, compute_residual):
         return newton_krylov(compute_residual, self._x,
                              f_tol=self.options['fatol'],
                              f_rtol=self.options['ftol'],
                              verbose=True, callback=self.callback)
 
 
 class DFSANESolver(ScipySolver):
     """
     Solve using a spectral residual jacobianless method
     """
     def __init__(self, residual, model, callback=None):
         ScipySolver.__init__(self, residual, model, callback)
-        self.options = {'ftol': 0, 'fatol': 1e-9}
 
     def scipy_solve(self, compute_residual):
         solution = root(compute_residual,
                         self._x,
                         method='df-sane',
                         options=self.options,
                         callback=self.callback)
         print("DF-SANE: {} ({} iterations, {})".format(
             solution.message,
             solution.nit,
             self.options))
         return solution.x.copy()
 
 
 class ToleranceManager:
     """
     Manage tolerance of non-linear solver
     """
     def __init__(self, start, end, rate, key='fatol'):
         self.start = start / rate  # just anticipating first multiplication
         self.end = end
         self.rate = rate
         self.key = key
 
     def __call__(self, cls):
         orig_init = cls.__init__
         orig_solve = cls.scipy_solve
         orig_updateState = cls.updateState
 
         def __init__(obj, *args, **kwargs):
             orig_init(obj, *args, **kwargs)
             obj.options[self.key] = self.start
 
         def scipy_solve(obj, *args, **kwargs):
             ftol = obj.options[self.key]
             ftol *= self.rate
 
             obj.options[self.key] = max(ftol, self.end)
             return orig_solve(obj, *args, **kwargs)
 
         def updateState(obj, *args, **kwargs):
             obj.options[self.key] = self.start
             return orig_updateState(obj, *args, **kwargs)
 
         cls.__init__ = __init__
         cls.scipy_solve = scipy_solve
         cls.updateState = updateState
         return cls