Fixing interpolation test

src/FiniteDifference/Derivative.hpp

@@ -12,6 +12,7 @@
 #define CENTRAL_B_ONE_SIDE 1
 #define FORWARD 2
 #define BACKWARD 3
+#define CENTRAL_SYM 4
 
 #include "util/mathutil.hpp"
 #include "Vector/map_vector.hpp"

src/FiniteDifference/FDScheme.hpp

@@ -234,7 +234,6 @@ private:
 			{
 				ke.eq = row - row_low * Sys_eqs::nvar;
 				ke.key = g_map.getGKey(it.get());
-				ke.key -= pd.getKP1();
 				return ke;
 			}
 
@@ -378,19 +377,22 @@ private:
 		if (g_dst.is_staggered() == true)
 			std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid must be normal " << std::endl;
 
-#ifdef SE_CLASS1
+		grid_key_dx<Grid_dst::dims> start;
+		grid_key_dx<Grid_dst::dims> stop;
 
-		if (g_map.getLocalDomainSize() != g_dst.getLocalDomainSize())
-			std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid in size does not match the mapping grid" << std::endl;
-#endif
+		for (size_t i = 0 ; i < Grid_dst::dims ; i++)
+		{
+			start.set_d(i,pd.getLow(i));
+			stop.set_d(i,g_map.size(i) - pd.getHigh(i));
+		}
 
 		// sub-grid iterator over the grid map
-		auto g_map_it = g_map.getDomainIterator();
+		auto g_map_it = g_map.getSubDomainIterator(start,stop);
 
 		// Iterator over the destination grid
 		auto g_dst_it = g_dst.getDomainIterator();
 
-		while (g_map_it.isNext() == true)
+		while (g_dst_it.isNext() == true)
 		{
 			typedef typename to_boost_vmpl<pos...>::type vid;
 			typedef boost::mpl::size<vid> v_size;
@@ -427,7 +429,88 @@ private:
 	 * \param it_d iterator that define where you want to impose
 	 *
 	 */
-	template<typename T, typename bop, typename iterator> void impose(const T & op ,
+	template<typename T, typename bop, typename iterator> void impose_dit(const T & op ,
+			                         bop num,
+									 long int id ,
+									 const iterator & it_d)
+	{
+		openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
+
+		auto it = it_d;
+		grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
+
+		std::unordered_map<long int,float> cols;
+
+		// iterate all the grid points
+		while (it.isNext())
+		{
+			// get the position
+			auto key = it.get();
+
+			// Add padding
+			for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
+				key.getKeyRef().set_d(i,key.getKeyRef().get(i) + pd.getLow(i));
+
+			// Calculate the non-zero colums
+			T::value(g_map,key,gs,spacing,cols,1.0);
+
+			// indicate if the diagonal has been set
+			bool is_diag = false;
+
+			// create the triplet
+			for ( auto it = cols.begin(); it != cols.end(); ++it )
+			{
+				trpl.add();
+				trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
+				trpl.last().col() = it->first;
+				trpl.last().value() = it->second;
+
+				if (trpl.last().row() == trpl.last().col())
+					is_diag = true;
+
+//				std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
+			}
+
+			// If does not have a diagonal entry put it to zero
+			if (is_diag == false)
+			{
+				trpl.add();
+				trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
+				trpl.last().col() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
+				trpl.last().value() = 0.0;
+			}
+
+			b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);
+
+			cols.clear();
+
+			// if SE_CLASS1 is defined check the position
+#ifdef SE_CLASS1
+//			T::position(key,gs,s_pos);
+#endif
+
+			++row;
+			++row_b;
+			++it;
+		}
+	}
+
+	/*! \brief Impose an operator
+	 *
+	 * This function impose an operator on a particular grid region to produce the system
+	 *
+	 * Ax = b
+	 *
+	 * ## Stokes equation 2D, lid driven cavity with one splipping wall
+	 * \snippet eq_unit_test.hpp Copy the solution to grid
+	 *
+	 * \param op Operator to impose (A term)
+	 * \param num right hand side of the term (b term)
+	 * \param id Equation id in the system that we are imposing
+	 * \param it_d iterator that define where you want to impose
+	 *
+	 */
+	template<typename T, typename bop, typename iterator> void impose_git(const T & op ,
 			                         bop num,
 									 long int id ,
 									 const iterator & it_d)
@@ -477,7 +560,6 @@ private:
 			b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);
 
 			cols.clear();
-//			std::cout << "\n";
 
 			// if SE_CLASS1 is defined check the position
 #ifdef SE_CLASS1
@@ -674,7 +756,7 @@ public:
 
         constant_b b(num);
 
-        impose(op,b,id,it);
+        impose_git(op,b,id,it);
 
 	}
 
@@ -693,14 +775,14 @@ public:
 	 * \param it_d iterator that define where you want to impose
 	 *
 	 */
-	template<typename T> void impose(const T & op ,
+	template<typename T> void impose_dit(const T & op ,
 									 typename Sys_eqs::stype num,
 									 long int id ,
 									 grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d)
 	{
 		constant_b b(num);
 
-		impose(op,b,id,it_d);
+		impose_dit(op,b,id,it_d);
 	}
 
 	/*! \brief Impose an operator
@@ -718,14 +800,14 @@ public:
 	 * \param it_d iterator that define where you want to impose
 	 *
 	 */
-	template<unsigned int prp, typename T, typename b_term, typename iterator> void impose(const T & op ,
+	template<unsigned int prp, typename T, typename b_term, typename iterator> void impose_dit(const T & op ,
 									 b_term & b_t,
 									 long int id ,
 									 const iterator & it_d)
 	{
 		grid_b<b_term,prp> b(b_t);
 
-		impose(op,b,id,it_d);
+		impose_dit(op,b,id,it_d);
 	}
 
 	//! type of the sparse matrix

src/FiniteDifference/Laplacian.hpp

@@ -125,4 +125,77 @@ public:
 };
 
 
+/*! \brief Laplacian second order approximation CENTRAL Scheme (with central derivative in the single)
+ *
+ * \verbatim
+
+          1.0
+           *
+           | -4.0
+   1.0 *---+---* 1.0
+           |
+           *
+          1.0
+
+ * \endverbatim
+ *
+ *
+ */
+template<typename arg, typename Sys_eqs>
+class Lap<arg,Sys_eqs,CENTRAL_SYM>
+{
+public:
+
+
+	/*! \brief Calculate which colums of the Matrix are non zero
+	 *
+	 * stub_or_real it is just for change the argument type when testing, in normal
+	 * conditions it is a distributed map
+	 *
+	 * \param pos position where the laplacian is calculated
+	 * \param gs Grid info
+	 * \param cols non-zero colums calculated by the function
+	 * \param coeff coefficent (constant in front of the derivative)
+	 *
+	 * ### Example
+	 *
+	 * \snippet FDScheme_unit_tests.hpp Laplacian usage
+	 *
+	 */
+	inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+	{
+		// for each dimension
+		for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
+		{
+			long int old_val = kmap.getKeyRef().get(i);
+			kmap.getKeyRef().set_d(i, kmap.getKeyRef().get(i) + 2);
+			arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[i]/spacing[i]/4.0);
+			kmap.getKeyRef().set_d(i,old_val);
+
+			old_val = kmap.getKeyRef().get(i);
+			kmap.getKeyRef().set_d(i, kmap.getKeyRef().get(i) - 2);
+			arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[i]/spacing[i]/4.0);
+			kmap.getKeyRef().set_d(i,old_val);
+
+			arg::value(g_map,kmap,gs,spacing,cols, - 2.0 * coeff/spacing[i]/spacing[i]/4.0);
+		}
+	}
+
+
+	/*! \brief Calculate the position where the derivative is calculated
+	 *
+	 * In case of non staggered case this function just return a null grid_key, in case of staggered,
+	 * the CENTRAL Laplacian scheme return the position of the staggered property
+	 *
+	 * \param position where we are calculating the derivative
+	 * \param gs Grid info
+	 * \param s_pos staggered position of the properties
+	 *
+	 */
+	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
+	{
+		return arg::position(pos,gs,s_pos);
+	}
+};
+
 #endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_LAPLACIAN_HPP_ */

src/Solvers/petsc_solver.hpp

@@ -696,8 +696,11 @@ class petsc_solver<double>
 
 public:
 
+	typedef Vector<double,PETSC_BASE> return_type;
+
 	~petsc_solver()
 	{
+		PETSC_SAFE_CALL(KSPDestroy(&ksp));
 	}
 
 	petsc_solver()
@@ -850,8 +853,14 @@ public:
 	 *
 	 *	\tparam impl Implementation of the SparseMatrix
 	 *
+	 * \param A sparse matrix
+	 * \param b vector
+	 * \param initial_guess true if x has the initial guess
+	 *
+	 * \return the solution
+	 *
 	 */
-	Vector<double,PETSC_BASE> solve(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & b)
+	Vector<double,PETSC_BASE> solve(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & b, bool initial_guess = false)
 	{
 		Mat & A_ = A.getMat();
 		const Vec & b_ = b.getVec();
@@ -864,6 +873,7 @@ public:
 
 		PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
 		PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
+		PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
 
 		Vector<double,PETSC_BASE> x(row,row_loc);
 		Vec & x_ = x.getVec();
@@ -876,6 +886,49 @@ public:
 		x.update();
 		return x;
 	}
+
+	/*! \brief Here we invert the matrix and solve the system
+	 *
+	 *  \warning umfpack is not a parallel solver, this function work only with one processor
+	 *
+	 *  \note if you want to use umfpack in a NON parallel, but on a distributed data, use solve with triplet
+	 *
+	 *	\tparam impl Implementation of the SparseMatrix
+	 *
+	 * \param A sparse matrix
+	 * \param b vector
+	 * \param x solution and initial guess
+	 * \param initial_guess true if x has the initial guess
+	 *
+	 * \return true if succeed
+	 *
+	 */
+	bool solve(SparseMatrix<double,int,PETSC_BASE> & A, Vector<double,PETSC_BASE> & x, const Vector<double,PETSC_BASE> & b)
+	{
+		Mat & A_ = A.getMat();
+		const Vec & b_ = b.getVec();
+
+		// We set the size of x according to the Matrix A
+		PetscInt row;
+		PetscInt col;
+		PetscInt row_loc;
+		PetscInt col_loc;
+
+		PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
+		PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
+		PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_TRUE));
+
+
+		Vec & x_ = x.getVec();
+
+		if (try_solve == true)
+			try_solve_simple(A_,b_,x_);
+		else
+			solve_simple(A_,b_,x_);
+
+		x.update();
+		return true;
+	}
 };
 
 #endif

src/Vector/Vector_petsc.hpp

@@ -100,7 +100,8 @@ class Vector<T,PETSC_BASE>
 	 */
 	void setPetsc() const
 	{
-		PETSC_SAFE_CALL(VecSetFromOptions(v));
+		if (v_created == false)
+		{PETSC_SAFE_CALL(VecSetType(v,VECMPI));}
 
 
 		// set the vector
@@ -366,6 +367,18 @@ public:
 
 		return *this;
 	}
+
+	/*! \brief Set to zero all the entries
+	 *
+	 *
+	 */
+	void setZero()
+	{
+		if (v_created == false)
+		{PETSC_SAFE_CALL(VecSetType(v,VECMPI));}
+
+		v_created = true;
+	}
 };
 
 

src/interpolation/interpolation.hpp

@@ -314,7 +314,7 @@ public:
 		}
 	}
 
-	template<unsigned int prp_v, unsigned int prp_g> void m2p(grid & gd, vector & vd)
+	template<unsigned int prp_g, unsigned int prp_v> void m2p(grid & gd, vector & vd)
 	{
 #ifdef SE_CLASS1
 

src/interpolation/interpolation_unit_tests.hpp

@@ -154,6 +154,10 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test )
 	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
 	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
 
+	// Do a ghost put
+
+	auto & v_cl = create_vcluster();
+
 	// We have to do a ghost get before interpolating m2p
 	// Before doing mesh to particle particle must be arranged
 	// into a grid like
@@ -191,14 +195,14 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test )
 	gd.ghost_get<0>();
 
 	grid_key_dx<2> start({3,3});
-	grid_key_dx<2> stop({(long int)gd.size(0) - 3,(long int)gd.size(1) - 3});
+	grid_key_dx<2> stop({(long int)gd.size(0) - 4,(long int)gd.size(1) - 4});
 
 	auto it6 = gd.getSubDomainIterator(start,stop);
 	while(it6.isNext())
 	{
 		auto key = it6.get();
 
-		gd.get<0>(key) = (double)rand()/RAND_MAX;;
+		gd.get<0>(key) = 5.0/*(double)rand()/RAND_MAX;*/;
 
 		++it6;
 	}
@@ -211,18 +215,36 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test )
 	momenta_grid_domain<decltype(gd),0>(gd,mg);
 	momenta_vector<decltype(vd),0>(vd,mv);
 
+	v_cl.sum(mg[0]);
+	v_cl.sum(mg[1]);
+	v_cl.sum(mv[0]);
+	v_cl.sum(mv[1]);
+	v_cl.execute();
+
 	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
 	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
 
 	momenta_grid_domain<decltype(gd),1>(gd,mg);
 	momenta_vector<decltype(vd),1>(vd,mv);
 
+	v_cl.sum(mg[0]);
+	v_cl.sum(mg[1]);
+	v_cl.sum(mv[0]);
+	v_cl.sum(mv[1]);
+	v_cl.execute();
+
 	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
 	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
 
 	momenta_grid_domain<decltype(gd),2>(gd,mg);
 	momenta_vector<decltype(vd),2>(vd,mv);
 
+	v_cl.sum(mg[0]);
+	v_cl.sum(mg[1]);
+	v_cl.sum(mv[0]);
+	v_cl.sum(mv[1]);
+	v_cl.execute();
+
 	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
 	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);