diff --git a/configure.ac b/configure.ac
index 51aaf839ee4d70f7fe299bdaf5bc0f9ef576e3a7..7675f169579547a69fce39b7f7f669e18afc66ea 100755
--- a/configure.ac
+++ b/configure.ac
@@ -136,8 +136,8 @@ AX_LIB_PETSC()
 AC_LANG_PUSH([C++])
 
 AC_CHECK_HEADER(quadmath.h, , [])
-AC_CHECK_LIB(quadmath, sinq, [ AC_CHECK_HEADER( quadmath.h ,[AC_DEFINE(HAVE_LIBQUADMATH,[],[Have quad math lib])
-				LIBQUADMATH=" -lquadmath "],[])
+AC_CHECK_LIB(quadmath, sinq, [  AC_DEFINE(HAVE_LIBQUADMATH,[],[Have quad math lib])
+				LIBQUADMATH=" -lquadmath "
 				], [])
 
 AC_LANG_POP([C++])
@@ -159,6 +159,14 @@ AX_LIB_HILBERT([],[echo "Cannot detect libhilbert, use the --with-libhilbert opt
 
 INCLUDES_PATH+=" -I/usr/local/include -I. -I../../openfpm_devices/src -I../../openfpm_data/src -I../../openfpm_io/src -I../../openfpm_vcluster/src -I../../src"
 
+####### Detect OpenMP
+
+AX_OPENMP([CFLAGS="$OPENMP_CFLAGS"
+           LDFLAGS="$OPENMP_LDFLAGS"],[])
+
+AC_SUBST(OPENMP_CFLAGS)
+AC_SUBST(OPENMP_LDFLAGS)
+
 ###### Check for se-class1
 
 AC_MSG_CHECKING(whether to build with security enhancement class1)
diff --git a/m4/ax_petsc_lib.m4 b/m4/ax_petsc_lib.m4
index ed14c831b4e230d19b4bd7c9b80b37fa7dd7a098..1bd159091052c10c40e9e9317a4cca0ec579c4ec 100755
--- a/m4/ax_petsc_lib.m4
+++ b/m4/ax_petsc_lib.m4
@@ -102,8 +102,10 @@ AC_DEFUN([AX_LIB_PETSC], [
                         old_CC=$CC
                         old_CFLAGS=$CFLAGS
                         old_LDFLAGS=$LDFLAGS
-                        CFLAGS="-I$with_petsc/include $HDF5_INCLUDE $METIS_INCLUDE "
-                        LDFLAGS="-L$with_petsc/lib $HDF5_LDFLAGS  $HDF5_LIBS $METIS_LIB -lmetis "
+			AX_OPENMP([CFLAGS="$OPENMP_CFLAGS"
+				   LDFLAGS="$OPENMP_LDFLAGS"],[])
+                        CFLAGS="$CFLAGS -I$with_petsc/include $HDF5_INCLUDE $METIS_INCLUDE "
+                        LDFLAGS="$LDFLAGS -L$with_petsc/lib $HDF5_LDFLAGS  $HDF5_LIBS $METIS_LIB -lmetis "
 			CC=$CXX
 
                         AC_LANG_SAVE
diff --git a/run.sh b/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..00fd49ba09b02f9470f470d6caafe2d33353fac6
--- /dev/null
+++ b/run.sh
@@ -0,0 +1,14 @@
+#! /bin/bash
+
+echo "RUN numerics test"
+
+source $HOME/openfpm_vars
+
+mpirun -np $3 ./src/numerics
+if [ $? -ne 0 ]; then
+   curl -X POST --data "payload={\"icon_emoji\": \":jenkins:\", \"username\": \"jenkins\"  , \"attachments\":[{ \"title\":\"Error:\", \"color\": \"#FF0000\", \"text\":\"$2 failed to complete the openfpm_numerics test \" }] }" https://hooks.slack.com/services/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
+   exit 1 ;
+fi
+
+curl -X POST --data "payload={\"icon_emoji\": \":jenkins:\", \"username\": \"jenkins\"  , \"attachments\":[{ \"title\":\"Info:\", \"color\": \"#00FF00\", \"text\":\"$2 completed succeffuly the openfpm_numerics test \" }] }" https://hooks.slack.com/services/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
+
diff --git a/src/FiniteDifference/Average.hpp b/src/FiniteDifference/Average.hpp
index 8e92b279e85a5f1071675bb1bae0e0f7da66667a..45867c6c7f0ea9ff238d9033c9d4af82eb747734 100644
--- a/src/FiniteDifference/Average.hpp
+++ b/src/FiniteDifference/Average.hpp
@@ -45,7 +45,7 @@ class Avg
 	 * In case of non staggered case this function just return a null grid_key, in case of staggered,
 	 *  it calculate how the operator shift the calculation in the cell
 	 *
-	 * \param position where we are calculating the derivative
+	 * \param pos position where we are calculating the derivative
 	 * \param gs Grid info
 	 * \param s_pos staggered position of the properties
 	 *
@@ -77,7 +77,7 @@ class Avg<d,arg,Sys_eqs,CENTRAL>
 	 * conditions it is a distributed map
 	 *
 	 * \param g_map It is the map explained in FDScheme
-	 * \param k_map position where the average is calculated
+	 * \param kmap position where the average is calculated
 	 * \param gs Grid info
 	 * \param cols non-zero colums calculated by the function
 	 * \param coeff coefficent (constant in front of the derivative)
@@ -113,7 +113,7 @@ class Avg<d,arg,Sys_eqs,CENTRAL>
 	 *
 	 * It follow the same concept of central derivative
 	 *
-	 * \param position where we are calculating the derivative
+	 * \param pos position where we are calculating the derivative
 	 * \param gs Grid info
 	 * \param s_pos staggered position of the properties
 	 *
@@ -160,7 +160,7 @@ class Avg<d,arg,Sys_eqs,FORWARD>
 	 * conditions it is a distributed map
 	 *
 	 * \param g_map It is the map explained in FDScheme
-	 * \param k_map position where the average is calculated
+	 * \param kmap position where the average is calculated
 	 * \param gs Grid info
 	 * \param cols non-zero colums calculated by the function
 	 * \param coeff coefficent (constant in front of the derivative)
@@ -188,7 +188,7 @@ class Avg<d,arg,Sys_eqs,FORWARD>
 	 * In case of non staggered case this function just return a null grid_key, in case of staggered,
 	 * the FORWARD scheme return the position of the staggered property
 	 *
-	 * \param position where we are calculating the derivative
+	 * \param pos position where we are calculating the derivative
 	 * \param gs Grid info
 	 * \param s_pos staggered position of the properties
 	 *
@@ -220,7 +220,7 @@ class Avg<d,arg,Sys_eqs,BACKWARD>
 	 * conditions it is a distributed map
 	 *
 	 * \param g_map It is the map explained in FDScheme
-	 * \param k_map position where the average is calculated
+	 * \param kmap position where the average is calculated
 	 * \param gs Grid info
 	 * \param cols non-zero colums calculated by the function
 	 * \param coeff coefficent (constant in front of the derivative)
@@ -247,7 +247,7 @@ class Avg<d,arg,Sys_eqs,BACKWARD>
 	 * In case of non staggered case this function just return a null grid_key, in case of staggered,
 	 * the BACKWARD scheme return the position of the staggered property
 	 *
-	 * \param position where we are calculating the derivative
+	 * \param pos position where we are calculating the derivative
 	 * \param gs Grid info
 	 * \param s_pos staggered position of the properties
 	 *
diff --git a/src/FiniteDifference/eq_unit_test_3d.hpp b/src/FiniteDifference/eq_unit_test_3d.hpp
index 83ec81d88f61477d81b41e7221efd270b42b1e86..f8b4bc3d03e10b21498bb0d94f74846d85858a75 100644
--- a/src/FiniteDifference/eq_unit_test_3d.hpp
+++ b/src/FiniteDifference/eq_unit_test_3d.hpp
@@ -22,57 +22,56 @@
 
 BOOST_AUTO_TEST_SUITE( eq_test_suite_3d )
 
-//!
-
+//! Specify the general caratteristic of system to solve
 struct lid_nn_3d_eigen
 {
-	// dimensionaly of the equation ( 3D problem ...)
+	//! dimensionaly of the equation ( 3D problem ...)
 	static const unsigned int dims = 3;
-	// number of fields in the system
+	//! number of fields in the system
 	static const unsigned int nvar = 4;
 
-	// boundary at X and Y
+	//! boundary at X and Y
 	static const bool boundary[];
 
-	// type of space float, double, ...
+	//! type of space float, double, ...
 	typedef float stype;
 
-	// type of base grid
+	//! type of base grid
 	typedef grid_dist_id<3,float,aggregate<float[3],float>,CartDecomposition<3,float>> b_grid;
 
-	// type of SparseMatrix for the linear solver
+	//! type of SparseMatrix for the linear solver
 	typedef SparseMatrix<double,int> SparseMatrix_type;
 
-	// type of Vector for the linear solver
+	//! type of Vector for the linear solver
 	typedef Vector<double> Vector_type;
 
-	// Define the underline grid is staggered
+	//! Define the underline grid is staggered
 	static const int grid_type = STAGGERED_GRID;
 };
 
 struct lid_nn_3d_petsc
 {
-	// dimensionaly of the equation ( 3D problem ...)
+	//! dimensionaly of the equation ( 3D problem ...)
 	static const unsigned int dims = 3;
-	// number of fields in the system
+	//! number of fields in the system
 	static const unsigned int nvar = 4;
 
-	// boundary at X and Y
+	//! boundary at X and Y
 	static const bool boundary[];
 
-	// type of space float, double, ...
+	//! type of space float, double, ...
 	typedef float stype;
 
-	// type of base grid
+	//! type of base grid
 	typedef grid_dist_id<3,float,aggregate<float[3],float>,CartDecomposition<3,float>> b_grid;
 
-	// type of SparseMatrix for the linear solver
+	//! type of SparseMatrix for the linear solver
 	typedef SparseMatrix<double,int,PETSC_BASE> SparseMatrix_type;
 
-	// type of Vector for the linear solver
+	//! type of Vector for the linear solver
 	typedef Vector<double,PETSC_BASE> Vector_type;
 
-	// Define the underline grid is staggered
+	//! Define the underline grid is staggered
 	static const int grid_type = STAGGERED_GRID;
 };
 
@@ -84,14 +83,13 @@ const bool lid_nn_3d_petsc::boundary[] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC
 // Constant Field
 struct eta
 {
+	//! define that eta is a constant field
 	typedef void const_field;
 
+	//! therutn the value of the constant
 	static float val()	{return 1.0;}
 };
 
-//#define SYSEQ_TYPE lid_nn_3d_eigen;
-//#include "Equations/stoke_flow_eq_3d.hpp"
-
 template<typename solver_type,typename lid_nn_3d> void lid_driven_cavity_3d()
 {
 	#include "Equations/stoke_flow_eq_3d.hpp"
diff --git a/src/Makefile.am b/src/Makefile.am
index 5b8f08cbb0171bbf37f466d111861903fbd74365..90ad2275f837baa87bb1a10154bd9869d37c2a94 100755
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,9 +1,9 @@
 
-LINKLIBS = $(LIBHILBERT_LIB) $(PETSC_LIB) $(SUITESPARSE_LIBS) $(LAPACK_LIBS) $(BLAS_LIBS)  $(METIS_LIB) $(PARMETIS_LIB) $(DEFAULT_LIB) $(PTHREAD_LIBS) $(OPT_LIBS) $(BOOST_LDFLAGS) $(BOOST_PROGRAM_OPTIONS_LIB) $(BOOST_IOSTREAMS_LIB) $(HDF5_LDFLAGS)  $(HDF5_LIBS) $(LIBQUADMATH)
+LINKLIBS = $(OPENMP_LDFLAGS) $(LIBHILBERT_LIB) $(PETSC_LIB) $(SUITESPARSE_LIBS) $(LAPACK_LIBS) $(BLAS_LIBS)  $(METIS_LIB) $(PARMETIS_LIB) $(DEFAULT_LIB) $(PTHREAD_LIBS) $(OPT_LIBS) $(BOOST_LDFLAGS) $(BOOST_PROGRAM_OPTIONS_LIB) $(BOOST_IOSTREAMS_LIB) $(HDF5_LDFLAGS)  $(HDF5_LIBS) $(LIBQUADMATH) $(OPENMP_LDFLAGS)
 
 noinst_PROGRAMS = numerics
 numerics_SOURCES = main.cpp ../../openfpm_vcluster/src/VCluster.cpp ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp
-numerics_CXXFLAGS = $(AM_CXXFLAGS) $(HDF5_CPPFLAGS) $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(PARMETIS_INCLUDE)  $(EIGEN_INCLUDE) -I/usr/local/libhilbert/include -Wno-deprecated-declarations -Wno-unused-local-typedefs
+numerics_CXXFLAGS = $(OPENMP_CFLAGS) $(LIBHILBERT_INCLUDE) $(AM_CXXFLAGS) $(HDF5_CPPFLAGS) $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(PARMETIS_INCLUDE) $(EIGEN_INCLUDE) $(PETSC_INCLUDE) -Wno-deprecated-declarations -Wno-unused-local-typedefs
 numerics_CFLAGS = $(CUDA_CFLAGS)
 numerics_LDADD = $(LINKLIBS) -lparmetis -lmetis
 nobase_include_HEADERS = Matrix/SparseMatrix.hpp Matrix/SparseMatrix_Eigen.hpp Matrix/SparseMatrix_petsc.hpp \
diff --git a/src/Matrix/SparseMatrix_petsc.hpp b/src/Matrix/SparseMatrix_petsc.hpp
index d129624d286ba4340ffa62f3cfc1a98c23416889..78d352d1ab1d7d4a3921231173b26324edec3724 100644
--- a/src/Matrix/SparseMatrix_petsc.hpp
+++ b/src/Matrix/SparseMatrix_petsc.hpp
@@ -24,22 +24,42 @@
 template<typename T>
 class triplet<T,PETSC_BASE>
 {
+	//! Row of the sparse matrix
 	PetscInt row_;
+
+	//! Colum of the sparse matrix
 	PetscInt col_;
+
+	//! Value of the Matrix
 	PetscScalar val_;
 
 public:
 
+	/*! \brief Return the row of the triplet
+	 *
+	 * \return the row index
+	 *
+	 */
 	PetscInt & row()
 	{
 		return row_;
 	}
 
+	/*! \brief Return the colum of the triplet
+	 *
+	 * \return the colum index
+	 *
+	 */
 	PetscInt & col()
 	{
 		return col_;
 	}
 
+	/*! \brief Return the value of the triplet
+	 *
+	 * \return the value
+	 *
+	 */
 	PetscScalar & value()
 	{
 		return val_;
diff --git a/src/Operators/Vector/vector_dist_operators.hpp b/src/Operators/Vector/vector_dist_operators.hpp
index 7adb79008e4f184796bd44a2120bdc3ffc018b1a..21ec3c8271f2e5dfb8f8f8f60df7162a18523715 100644
--- a/src/Operators/Vector/vector_dist_operators.hpp
+++ b/src/Operators/Vector/vector_dist_operators.hpp
@@ -499,12 +499,15 @@ public:
 template<unsigned int prp>
 class vector_dist_expression<prp,float>
 {
+	//! constant value
 	float d;
 
 public:
 
+	//! type of object the structure return then evaluated
 	typedef float vtype;
 
+	//! constrictir from constant value
 	inline vector_dist_expression(const float & d)
 	:d(d)
 	{}
@@ -519,7 +522,9 @@ public:
 
 	/*! \brief Evaluate the expression
 	 *
-	 * It just return the velue set in the constructor
+	 * It just return the value set in the constructor
+	 *
+	 * \return the constant value set in the constructor
 	 *
 	 */
 	inline float value(const vect_dist_key_dx & k) const
diff --git a/src/Operators/Vector/vector_dist_operators_functions.hpp b/src/Operators/Vector/vector_dist_operators_functions.hpp
index c0d9dfb3fd9079be6597c13515ca397619af21c9..ff7308e14a8ded59da50a283a9d731fb18e8edd2 100644
--- a/src/Operators/Vector/vector_dist_operators_functions.hpp
+++ b/src/Operators/Vector/vector_dist_operators_functions.hpp
@@ -49,14 +49,6 @@ fun_name(const vector_dist_expression_op<exp1,exp2_,op1> & va)\
 	return exp_sum;\
 }\
 \
-template<typename exp1, typename exp2_, unsigned int op1>\
-inline vector_dist_expression_op<vector_dist_expression<0,double>,void,OP_ID>\
-fun_name(double d)\
-{\
-	vector_dist_expression_op<vector_dist_expression<0,double>,void,OP_ID> exp_sum( (vector_dist_expression<0,double>(d)) );\
-\
-	return exp_sum;\
-}\
 \
 template<unsigned int prp1, typename v1>\
 inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,void,OP_ID>\