Added stencil grid + documentation

CHANGELOG.md

@@ -14,6 +14,8 @@ All notable changes to this project will be documented in this file.
 - Raw reader for grid (see ...)
 - A way to specify names for proeprties and select properties to write
 - Ghost put on grid
+- getDomainIterator stencil for faster stencil codes iterators
+- Agebraic multigrid solvers interface for linear systems
 
 ### Fixed
 - Installation of PETSC in case with MUMPS try without MUMPS
@@ -22,9 +24,11 @@ All notable changes to this project will be documented in this file.
 - Bug in VTK writer binary: long int are not supported removing output
 - Bug in FDScheme in the constructor with stencil bigger than one
 - Bug Fixed Memory leak in petsc solver
+- Bug Performance bug in the grid iterator
 
 ### Changed
-- CellList types has changed
+- CellList types has changed (one additional template parameter)
+- Gris iterator types has changes (one additional template parameter)
 
 ## [0.8.0] 28 February 2016
 

example/Grid/3_gray_scott_3d/main.cpp

@@ -140,6 +140,9 @@ int main(int argc, char* argv[])
 	double uFactor = deltaT * du/(spacing[x]*spacing[x]);
 	double vFactor = deltaT * dv/(spacing[x]*spacing[x]);
 
+	timer tot_sim;
+	tot_sim.start();
+
 	for (size_t i = 0; i < timeSteps; ++i)
 	{
 		if (i % 300 == 0)
@@ -199,6 +202,9 @@ int main(int argc, char* argv[])
 		}
 	}
 	
+	tot_sim.stop();
+	std::cout << "Total simulation: " << tot_sim.getwct() << std::endl;
+
 	//! \cond [time stepping] \endcond
 
 	/*!

example/Numerics/Vortex_in_cell/main_vic_petsc.cpp

@@ -116,19 +116,19 @@ typedef grid_dist_id<3,float,aggregate<float[3]>> grid_type;
 typedef vector_dist<3,float,aggregate<float[3],float[3],float[3],float[3],float[3]>> particles_type;
 
 // radius of the torus
-float ringr1 = 5.0/4.0;
+float ringr1 = 1.0;
 // radius of the core of the torus
-float ringr2 = 0.5;
+float sigma = 1.0/3.523;
 // Reynold number
-float tgtre  = 10000.0;
+float tgtre  = 7500.0;
 // Noise factor for the ring vorticity on z
-float ringnz = 0.00009;
+float ringnz = 0.01;
 
 // Kinematic viscosity
-float nu = 0.0001535;
+float nu = 1.0/tgtre;
 
 // Time step
-float dt = 0.006;
+float dt = 0.025;
 
 // All the properties by index
 constexpr unsigned int vorticity = 0;
@@ -225,13 +225,13 @@ void init_ring(grid_type & gr, const Box<3,float> & domain)
 
 	for (size_t i = 0 ; i < nk ; i++)
 	{
-	     ak[i] = 0.0/*rand()/RAND_MAX*/;
-	     bk[i] = 0.0/*rand()/RAND_MAX*/;
+	     ak[i] = rand()/RAND_MAX;
+	     bk[i] = rand()/RAND_MAX;
 	}
 
 	// We calculate the circuitation gamma
 	float gamma = nu * tgtre;
-	float rinv2 = 1.0f/(ringr2*ringr2);
+	float rinv2 = 1.0f/(sigma*sigma);
 	float max_vorticity = gamma*rinv2/M_PI;
 
 	// We go through the grid to initialize the vortex
@@ -247,27 +247,20 @@ void init_ring(grid_type & gr, const Box<3,float> & domain)
         float tz = (key.get(z)-2)*gr.spacing(z) + domain.getLow(z);
         float theta1 = atan2((ty-2.5f),(tz-2.5f));
 
+
+
         float noise = 0.0f;
-        for (int kk=1 ; kk < nk; kk++)
-        	noise = noise + sin(kk*(theta1+2.0f*M_PI*ak[kk])) + cos(kk*(theta1+2.0f*M_PI*bk[kk]));
+ //       for (int kk=1 ; kk < nk; kk++)
+ //       	noise = noise + sin(kk*(theta1+2.0f*M_PI*ak[kk])) + cos(kk*(theta1+2.0f*M_PI*bk[kk]));
 
         float rad1r  = sqrt((ty-2.5f)*(ty-2.5f) + (tz-2.5f)*(tz-2.5f)) - ringr1*(1.0f + ringnz * noise);
-        float rad1t = tx - 2.5f;
+        float rad1t = tx - 1.0f;
         float rad1sq = rad1r*rad1r + rad1t*rad1t;
         float radstr = -exp(-rad1sq*rinv2)*rinv2*gamma/M_PI;
         gr.template get<vorticity>(key_d)[x] = 0.0f;
         gr.template get<vorticity>(key_d)[y] = -radstr * cos(theta1);
         gr.template get<vorticity>(key_d)[z] = radstr * sin(theta1);
 
-/*        theta1 = atan2((ty-2.5f),(tz-2.5f));
-        rad1r  = sqrt((ty-2.5f)*(ty-2.5f) + (tz-2.5f)*(tz-2.5f)) + ringr1*(1.0f + ringnz * noise);
-        rad1t = tx - 2.5f;
-        rad1sq = rad1r*rad1r + rad1t*rad1t;
-        float rad1sqTILDA = rad1sq*rinv2;
-        radstr = exp(-rad1sqTILDA)*rinv2*gamma/M_PI;
-        gr.template get<vorticity>(key_d)[x] = 0.0f;
-        gr.template get<vorticity>(key_d)[y] = gr.template get<vorticity>(key_d)[y] + radstr * cos(theta1);
-        gr.template get<vorticity>(key_d)[z] = gr.template get<vorticity>(key_d)[z] - radstr * sin(theta1);*/
 
 		++it;
 	}
@@ -494,8 +487,6 @@ void helmotz_hodge_projection(grid_type & gr, const Box<3,float> & domain)
 	// Set the maximum number of iterations
 	solver.setMaxIter(500);
 
-	solver.log_monitor();
-
 	timer tm_solve;
 	tm_solve.start();
 
@@ -704,7 +695,6 @@ void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc
 		solver.setSolver(KSPBCGS);
 		solver.setAbsTol(0.01);
 		solver.setMaxIter(500);
-		solver.log_monitor();
 
 		// Get the sparse matrix that represent the left-hand-side
 		// of the equation
@@ -1153,7 +1143,7 @@ int main(int argc, char* argv[])
 	openfpm_init(&argc,&argv);
 
 	// Domain, a rectangle
-	Box<3,float> domain({0.0,0.0,0.0},{22.0,5.0,5.0});
+	Box<3,float> domain({0.0,0.0,0.0},{22.0,5.57,5.57});
 
 	// Ghost (Not important in this case but required)
 	Ghost<3,long int> g(2);
@@ -1239,7 +1229,7 @@ int main(int argc, char* argv[])
 	}
 
 	// Time Integration
-	for ( ; i < 10 ; i++)
+	for ( ; i < 10001 ; i++)
 	{
 		// do step 4-5-6-7
 		do_step(particles,g_vort,g_vel,g_dvort,domain,inte,phi_s);
@@ -1258,6 +1248,9 @@ int main(int argc, char* argv[])
 		inte.template p2m<vorticity,vorticity>(particles,g_vort);
 		g_vort.template ghost_put<add_,vorticity>();
 
+		// helmotz-hodge projection
+		helmotz_hodge_projection(g_vort,domain);
+
 		remesh(particles,g_vort,domain);
 
 		// print the step number

openfpm_data @ 8ac02bf2

-Subproject commit dfaf8538f94caf51229ec40bc51af23ad8bea934
+Subproject commit 8ac02bf29a070237ca028a9ecddbd7ed81a4457d

openfpm_io @ bcc455d1

-Subproject commit 223f739a00a1d92034e71b1aa821f244c68f4ced
+Subproject commit bcc455d1dd01285b2f5f78497566720e7d8c0bf5

openfpm_numerics @ 58c0aea7

-Subproject commit 43a521cb99e7859af4392e70e0eb6196646f31d4
+Subproject commit 58c0aea7c21597e43967f700b54db58fe90b941e

openfpm_pdata.doc

@@ -1895,7 +1895,7 @@ ENABLE_PREPROCESSING   = YES
 # The default value is: NO.
 # This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
 
-MACRO_EXPANSION        = NO
+MACRO_EXPANSION        = YES
 
 # If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES then
 # the macro expansion is limited to the macros specified with the PREDEFINED and
@@ -1903,7 +1903,7 @@ MACRO_EXPANSION        = NO
 # The default value is: NO.
 # This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
 
-EXPAND_ONLY_PREDEF     = NO
+EXPAND_ONLY_PREDEF     = YES
 
 # If the SEARCH_INCLUDES tag is set to YES the includes files in the
 # INCLUDE_PATH will be searched if a #include is found.
@@ -1935,7 +1935,7 @@ INCLUDE_FILE_PATTERNS  =
 # recursively expanded use the := operator instead of the = operator.
 # This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
 
-PREDEFINED             =
+PREDEFINED             = HAVE_PETSC
 
 # If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then this
 # tag can be used to specify a list of macro names that should be expanded. The

openfpm_vcluster @ 09caef4a

-Subproject commit 7030e7f76e5d80a297e4cc9b2d4f3ad98beb8d0f
+Subproject commit 09caef4afa48f5b4013328562b388821be422993

src/Decomposition/CartDecomposition.hpp

@@ -1388,12 +1388,9 @@ public:
 
 	/*! \brief Get the domain Cells
 	 *
-	 * It perform a linearization of the domain cells using the extension provided in gs
+	 * It return all the cells-id that are inside the processor-domain
 	 *
-	 *
-	 * \param shift Cell padding
-	 * \param cell_shift where the domain cell start
-	 * \param gs grid extension
+	 * \return the cells id inside the domain
 	 *
 	 */
 	openfpm::vector<size_t> & getDomainCells()
@@ -1401,15 +1398,26 @@ public:
 		return domain_nn_calculator_cart<dim>::getDomainCells();
 	}
 
-	/*! \brief Get the CRS domain Cells
+	/*! \brief Get the CRS domain Cells with normal neighborhood
+	 *
+	 * In case of symmetric interaction the neighborhood cells of
+	 * a cell is different
+	 *
+	 * \verbatim
+
+	   Symmetric      Normal
+
+	    * * *         * * *
+	      X *         * X *
+	                  * * *
+	  \endverbatim
 	 *
-	 * It perform a linearization of the domain cells using the extension provided in gs
 	 *
+	 * In case of CRS scheme some cells has the symmetric neighborhood
+	 * some others has more complex neighborhood. This function return
+	 * all the cells with normal neighborhood
 	 *
-	 * \param shift Cell padding
-	 * \param cell_shift where the domain cell start
-	 * \param gs grid extension
-	 * \param loc_box processor sub-domains
+	 * \return the cell-id of the cells inside the processor-domain with normal neighborhood
 	 *
 	 */
 	openfpm::vector<size_t> & getCRSDomainCells()
@@ -1430,13 +1438,8 @@ public:
 
 	/*! \brief Get the CRS anomalous cells
 	 *
-	 * This function include also a linearization of the indexes
+	 * This function return the anomalous cells
 	 *
-	 * \param shift Shifting point
-	 * \param cell_shift where the processor cell-list start (In case of symmetric)
-	 *                   exist one global cell-list, but each processor span only one
-	 *                   part of it
-	 * \param gs grid extension
 	 *
 	 * \return the anomalous cells with neighborhood
 	 *

src/Decomposition/Distribution/parmetis_util.hpp

@@ -506,7 +506,9 @@ public:
 
 	/*! \brief Copy the object
 	 *
-	 * \param object to copy
+	 * \param pm object to copy
+	 *
+	 * \return itself
 	 *
 	 */
 	const Parmetis<Graph> & operator=(const Parmetis<Graph> & pm)
@@ -523,7 +525,9 @@ public:
 
 	/*! \brief Copy the object
 	 *
-	 * \param object to copy
+	 * \param pm object to copy
+	 *
+	 * \return itself
 	 *
 	 */
 	const Parmetis<Graph> & operator=(Parmetis<Graph> && pm)

src/dec_optimizer.hpp → src/Decomposition/dec_optimizer.hpp

@@ -144,7 +144,7 @@ private:
 	template<unsigned int p_sub> void fill_domain(Graph & graph,const Box<dim,size_t> & box, long int ids)
 	{
 		// Create a subgrid iterator
-		grid_key_dx_iterator_sub<dim,do_not_print_warning_on_adjustment<dim>> g_sub(gh,box.getKP1(),box.getKP2());
+		grid_key_dx_iterator_sub<dim,no_stencil,do_not_print_warning_on_adjustment<dim>> g_sub(gh,box.getKP1(),box.getKP2());
 
 		// iterate through all grid points
 
@@ -207,7 +207,7 @@ private:
 		for (size_t d = 0 ; d < v_w.size() ; d++)
 		{
 			// Create a sub-grid iterator
-			grid_key_dx_iterator_sub_bc<dim,do_not_print_warning_on_adjustment<dim>> g_sub(gh,v_w.template get<wavefront<dim>::start>(d),v_w.template get<wavefront<dim>::stop>(d),bc);
+			grid_key_dx_iterator_sub_bc<dim,no_stencil,do_not_print_warning_on_adjustment<dim>> g_sub(gh,v_w.template get<wavefront<dim>::start>(d),v_w.template get<wavefront<dim>::stop>(d),bc);
 
 			// iterate through all grid points
 
@@ -296,7 +296,7 @@ private:
 				// Create an iterator of the expanded wavefront
 				grid_key_dx<dim> start = grid_key_dx<dim>(v_w.template get<wavefront<dim>::start>(d)) + w_comb[d];
 				grid_key_dx<dim> stop = grid_key_dx<dim>(v_w.template get<wavefront<dim>::stop>(d)) + w_comb[d];
-				grid_key_dx_iterator_sub<dim,do_not_print_warning_on_adjustment<dim>> it(gh,start,stop);
+				grid_key_dx_iterator_sub<dim,no_stencil,do_not_print_warning_on_adjustment<dim>> it(gh,start,stop);
 
 				// for each sub-domain in the expanded wavefront
 				while (it.isNext())

src/Grid/Iterators/grid_dist_id_iterator.hpp

@@ -25,7 +25,7 @@
  * \tparam impl implementation
  *
  */
-template<unsigned int dim, typename device_grid, int impl >
+template<unsigned int dim, typename device_grid, int impl,typename stencil = no_stencil >
 class grid_dist_iterator
 {
 
@@ -41,8 +41,8 @@ class grid_dist_iterator
  * \tparam impl implementation
  *
  */
-template<unsigned int dim, typename device_grid>
-class grid_dist_iterator<dim,device_grid,FREE>
+template<unsigned int dim, typename device_grid, typename stencil>
+class grid_dist_iterator<dim,device_grid,FREE,stencil>
 {
 	//! grid list counter
 	size_t g_c;
@@ -54,7 +54,7 @@ class grid_dist_iterator<dim,device_grid,FREE>
 	const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext;
 
 	//! Actual iterator
-	grid_key_dx_iterator_sub<dim> a_it;
+	grid_key_dx_iterator_sub<dim,stencil> a_it;
 
 	//! stop point (is the grid size)
 	grid_key_dx<dim> stop;
@@ -91,6 +91,26 @@ class grid_dist_iterator<dim,device_grid,FREE>
 		selectValidGrid();
 	}
 
+	/*! \brief Constructor of the distributed grid iterator with
+	 *         stencil support
+	 *
+	 * \param gk std::vector of the local grid
+	 * \param gdb_ext set of local subdomains
+	 * \param stop end point
+	 * \param stencil_pnt stencil points
+	 *
+	 */
+	grid_dist_iterator(const openfpm::vector<device_grid> & gk,
+			           const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext,
+					   const grid_key_dx<dim> & stop,
+					   const grid_key_dx<dim> (& stencil_pnt)[stencil::nsp])
+	:g_c(0),gList(gk),gdb_ext(gdb_ext),a_it(stencil_pnt),stop(stop)
+	{
+		// Initialize the current iterator
+		// with the first grid
+		selectValidGrid();
+	}
+
 	// Destructor
 	~grid_dist_iterator()
 	{
@@ -102,7 +122,7 @@ class grid_dist_iterator<dim,device_grid,FREE>
 	 *
 	 */
 
-	inline grid_dist_iterator<dim,device_grid,FREE> operator++()
+	inline grid_dist_iterator<dim,device_grid,FREE,stencil> & operator++()
 	{
 		++a_it;
 
@@ -191,6 +211,8 @@ class grid_dist_iterator<dim,device_grid,FREE>
 	 * \see grid_dist_key_dx
 	 * \see grid_dist_iterator
 	 *
+	 * \param k key position in local coordinates
+	 *
 	 * \return the global position in the grid
 	 *
 	 */
@@ -206,6 +228,18 @@ class grid_dist_iterator<dim,device_grid,FREE>
 
 		return k_glob;
 	}
+
+	/*! \brief Return the stencil point offset
+	 *
+	 * \tparam id
+	 *
+	 * \return linearized distributed key
+	 *
+	 */
+	template<unsigned int id> inline grid_dist_lin_dx getStencil()
+	{
+		return grid_dist_lin_dx(g_c,a_it.getStencil<id>());
+	}
 };
 
 
@@ -218,8 +252,8 @@ class grid_dist_iterator<dim,device_grid,FREE>
  * \tparam impl implementation
  *
  */
-template<unsigned int dim, typename device_grid>
-class grid_dist_iterator<dim,device_grid,FIXED>
+template<unsigned int dim, typename device_grid,typename stencil>
+class grid_dist_iterator<dim,device_grid,FIXED,stencil>
 {
 	//! grid list counter
 	size_t g_c;
@@ -231,7 +265,7 @@ class grid_dist_iterator<dim,device_grid,FIXED>
 	const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext;
 
 	//! Actual iterator
-	grid_key_dx_iterator<dim> a_it;
+	grid_key_dx_iterator<dim,stencil> a_it;
 
 	/*! \brief from g_c increment g_c until you find a valid grid
 	 *
@@ -289,7 +323,7 @@ class grid_dist_iterator<dim,device_grid,FIXED>
 	 *
 	 */
 
-	grid_dist_iterator<dim,device_grid,FIXED> operator++()
+	grid_dist_iterator<dim,device_grid,FIXED> &  operator++()
 	{
 		++a_it;
 
@@ -349,6 +383,8 @@ class grid_dist_iterator<dim,device_grid,FIXED>
 	 * \see grid_dist_key_dx
 	 * \see grid_dist_iterator
 	 *
+	 * \param k local coordinates to convert into global
+	 *
 	 * \return the global position in the grid
 	 *
 	 */
@@ -364,6 +400,18 @@ class grid_dist_iterator<dim,device_grid,FIXED>
 
 		return k_glob;
 	}
+
+	/*! \brief Return the stencil point offset
+	 *
+	 * \tparam id
+	 *
+	 * \return linearized distributed key
+	 *
+	 */
+	template<unsigned int id> inline grid_dist_lin_dx getStencil()
+	{
+		return grid_dist_lin_dx(g_c,a_it.getStencil<id>());
+	}
 };
 
 #endif /* GRID_DIST_ID_ITERATOR_SUB_HPP_ */

src/Grid/Iterators/grid_dist_id_iterator_dec.hpp

@@ -161,7 +161,7 @@ class grid_dist_id_iterator_dec
 	 *
 	 */
 
-	inline grid_dist_id_iterator_dec<Decomposition> operator++()
+	inline grid_dist_id_iterator_dec<Decomposition> & operator++()
 	{
 		++a_it;
 

src/Grid/Iterators/grid_dist_id_iterator_sub.hpp

@@ -26,7 +26,7 @@ class grid_dist_iterator_sub
 	{
 		//! start point where iterate
 		grid_key_dx<dim> start;
-		// ! stop point where iterate
+		//! stop point where iterate
 		grid_key_dx<dim> stop;
 	};
 
@@ -165,7 +165,7 @@ class grid_dist_iterator_sub
 	 *
 	 */
 
-	inline grid_dist_iterator_sub<dim,device_grid> operator++()
+	inline grid_dist_iterator_sub<dim,device_grid> & operator++()
 	{
 		++a_it;
 

src/Grid/Iterators/grid_dist_id_iterators_unit_tests.hpp

@@ -281,6 +281,98 @@ void Test3D_decit(const Box<3,float> & domain, long int k)
 	}
 }
 
+void Test3D_stencil(const Box<3,float> & domain, long int k)
+{
+	grid_key_dx<3> star_stencil_3D[7] = {{0,0,0},
+	                                         {0,0,-1},
+											 {0,0,1},
+											 {0,-1,0},
+											 {0,1,0},
+											 {-1,0,0},
+											 {1,0,0}};
+
+	{
+		Vcluster & v_cl = create_vcluster();
+
+		if ( v_cl.getProcessingUnits() > 32 )
+			return;
+
+		long int big_step = k / 30;
+		big_step = (big_step == 0)?1:big_step;
+		long int small_step = 21;
+
+		print_test( "Testing grid stencil iterator k<=",k);
+
+		// 3D test
+		for ( ; k >= 2 ; k-= (k > 2*big_step)?big_step:small_step )
+		{
+			BOOST_TEST_CHECKPOINT( "Testing grid skin iterator from decomposition k<=" << k );
+
+			// grid size
+			size_t sz[3];
+			sz[0] = k;
+			sz[1] = k;
+			sz[2] = k;
+
+			if (k <= 9)
+				continue;
+
+			Ghost<3,long int> g(1);
+
+			// Distributed grid with id decomposition
+			grid_dist_id<3, float, aggregate<long int>, CartDecomposition<3,float>> g_dist(sz,domain,g);
+
+			// fill the grid with values
+
+			auto it = g_dist.getDomainGhostIterator();
+
+			while (it.isNext())
+			{
+				auto p = it.get();
+				auto gkey = it.getGKey(p);
+
+				g_dist.template get<0>(p) = gkey.get(0) + gkey.get(1) + gkey.get(2);
+
+				++it;
+			}
+
+			g_dist.ghost_get<0>();
+
+			auto st_it = g_dist.getDomainIteratorStencil(star_stencil_3D);
+
+			bool ret = true;
+
+			while (st_it.isNext())
+			{
+				// center point
+				auto Cp = st_it.getStencil<0>();
+
+				// plus,minus X,Y,Z
+				auto mx = st_it.getStencil<1>();
+				auto px = st_it.getStencil<2>();
+				auto my = st_it.getStencil<3>();
+				auto py = st_it.getStencil<4>();
+				auto mz = st_it.getStencil<5>();
+				auto pz = st_it.getStencil<6>();
+
+				size_t sum = 6*g_dist.template get<0>(Cp) -
+						     g_dist.template get<0>(mx) -
+							 g_dist.template get<0>(px) -
+							 g_dist.template get<0>(my) -
+							 g_dist.template get<0>(py) -
+							 g_dist.template get<0>(mz) -
+							 g_dist.template get<0>(pz);
+
+				ret &= (sum == 0);
+
+				++st_it;
+			}
+
+			BOOST_REQUIRE_EQUAL(ret,true);
+		}
+
+	}
+}
 
 // Test decomposition grid iterator
 
@@ -396,6 +488,17 @@ BOOST_AUTO_TEST_CASE( grid_dist_id_decomposition_iterator )
 	Test3D_decit(domain3,k);
 }
 
+BOOST_AUTO_TEST_CASE( grid_dist_id_iterator_stencil )
+{
+	// Domain
+	Box<3,float> domain3({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+	size_t k = 128*128*128*create_vcluster().getProcessingUnits();
+	k = std::pow(k, 1/3.);
+	Test3D_stencil(domain3,k);
+}
+
+
 BOOST_AUTO_TEST_CASE( grid_dist_it_iterators_skin_test )
 {