Last Finite Differences not working required several other components

src/FiniteDifference/Derivative.hpp

@@ -10,6 +10,7 @@
 
 #define CENTRAL 0
 #define CENTRAL_B_ONE_SIDE 1
+#define FORWARD 2
 
 #include "util/mathutil.hpp"
 #include "Vector/map_vector.hpp"
@@ -63,6 +64,13 @@ class D<d,arg,Sys_eqs,CENTRAL>
 	 */
 	inline static void value(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
 	{
+		// if the system is staggered the CENTRAL derivative is equivalent to a forward derivative
+		if (is_grid_staggered<Sys_eqs>::value() == true)
+		{
+			D<d,arg,Sys_eqs,FORWARD>::value(pos,gs,cols,coeff);
+			return;
+		}
+
 		// forward
 		if (Sys_eqs::boundary[d] == PERIODIC )
 		{
@@ -227,5 +235,57 @@ public:
 };
 
 
+/*! \brief Derivative FORWARD on direction i
+ *
+ *
+ */
+template<unsigned int d, typename arg, typename Sys_eqs>
+class D<d,arg,Sys_eqs,FORWARD>
+{
+	public:
+
+	/*! \brief fill the row
+	 *
+	 *
+	 */
+	inline static void value(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
+	{
+		// forward
+		if (Sys_eqs::boundary[d] == PERIODIC )
+		{
+			long int old_val = pos.get(d);
+			pos.set_d(d, openfpm::math::positive_modulo(pos.get(d) + 1, gs.size(d)));
+			arg::value(pos,gs,cols,coeff);
+			pos.set_d(d,old_val);
+		}
+		else
+		{
+			long int old_val = pos.get(d);
+			pos.set_d(d, pos.get(d) + 1);
+			arg::value(pos,gs,cols,coeff);
+			pos.set_d(d,old_val);
+		}
+
+		// backward
+		arg::value(pos,gs,cols,-coeff);
+	}
+
+
+	/*! \brief Calculate the position where the derivative is calculated
+	 *
+	 * In case on non staggered case this function just return pos, in case of staggered,
+	 *  it calculate where the operator is calculated on a staggered grid
+	 *
+	 *  \param pos from the position
+	 *  \param fld Field we are deriving, if not provided the function just return pos
+	 *  \param s_pos position of the properties in the staggered grid
+	 *
+	 */
+	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, long int fld = -1, const openfpm::vector<comb<Sys_eqs::dims>> & s_pos = openfpm::vector<comb<Sys_eqs::dims>>())
+	{
+		return pos;
+	}
+};
+
 
 #endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_ */

src/FiniteDifference/FDScheme.hpp

@@ -22,19 +22,106 @@
 template<typename Sys_eqs>
 class FDScheme
 {
+	// Padding
+	Padding<Sys_eqs::dims> pd;
+
+	// Sparse Matrix
 	openfpm::vector<triplet<typename Sys_eqs::stype>> trpl;
 
+	openfpm::vector<typename Sys_eqs::stype> b;
+
+	// Domain Grid informations
+	const grid_sm<Sys_eqs::dims,void> & gs;
+
+	// mapping grid
+	grid_dist_id<Sys_eqs::dims,Sys_eqs::stype,scalar<size_t>,Sys_eqs::grid_type::decomposition> g_map;
+
+	// row of the matrix
+	size_t row;
+
+	// row on b
+	size_t row_b;
+
+	/*! \brief Check if the Matrix is consistent
+	 *
+	 */
+	void consistency()
+	{
+		// A and B must have the same rows
+		if (row != row_b)
+			std::cerr << "Error " << __FILE__ << ":" << __LINE__ << "the term B and the Matrix A for Ax=B must contain the same number of rows\n";
+
+		// Indicate all the non zero rows
+		openfpm::vector<bool> nz_rows;
+		nz_rows.resize(row_b);
+
+		for (size_t i = 0 ; i < trpl.size() ; i++)
+		{
+			nz_rows.get(trpl.get(i).i) = true;
+		}
+
+		// Indicate all the non zero colums
+		openfpm::vector<bool> nz_cols;
+		for (size_t i = 0 ; i < trpl.size() ; i++)
+		{
+			nz_cols.get(trpl.get(i).j) = true;
+		}
+
+		// all the rows must have a non zero element
+		for (size_t i = 0 ; i < nz_rows.size() ; i++)
+		{
+			if (nz_rows.get(i) == false)
+				std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix not all the row are filled\n";
+		}
+
+		// all the colums must have a non zero element
+		for (size_t i = 0 ; i < nz_cols.size() ; i++)
+		{
+			if (nz_cols.get(i) == false)
+				std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix not all the row are filled\n";
+		}
+	}
+
 public:
 
+	/*! \brief Constructor
+	 *
+	 * \param pd Padding
+	 * \param gs grid infos where Finite differences work
+	 *
+	 */
+	FDScheme(Padding<Sys_eqs::dims> & pd, const grid_sm<Sys_eqs::dims,void> & gs)
+	:pd(pd),gs(gs)
+	{
+		// Counter
+		size_t cnt = 0;
+
+		// Resize the b term
+		b.resize(gs.size());
+
+		// Create the re-mapping-grid
+		auto it = g_map.getDomainIterator();
+
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			g_map.get(key) = cnt;
+
+			++cnt;
+			++it;
+		}
+
+
+	}
+
 	/*! \brief Impose an operator
 	 *
 	 *
 	 *
 	 */
-	template<typename T> void impose(T & op, const grid_sm<Sys_eqs::dims,void> & gs , grid_dist_iterator_sub<Sys_eqs::dims,typename Sys_eqs::b_grid::d_grid> it)
+	template<typename T> void imposeA(const T & op , grid_dist_iterator_sub<Sys_eqs::dims,typename Sys_eqs::b_grid::d_grid> it)
 	{
-		size_t lin = 0;
-
 		std::unordered_map<long int,float> cols;
 
 		// iterate all the grid points
@@ -54,7 +141,7 @@ public:
 			for ( auto it = cols.begin(); it != cols.end(); ++it )
 			{
 				trpl.add();
-				trpl.last().i = lin;
+				trpl.last().i = row;
 				trpl.last().j = it->first;
 				trpl.last().value = it->second;
 
@@ -64,25 +151,29 @@ public:
 			cols.clear();
 			std::cout << "\n";
 
+			++row;
 			++it;
 		}
 	}
 
-	/*! \brief produce the Matrix
+	/*! \brief Impose an operator
+	 *
 	 *
-	 *  \tparam Syst_eq System of equations, or a single equation
 	 *
 	 */
-	template<typename Grid> SparseMatrix<typename Sys_eqs::stype> getMatrix(const Grid & g)
+	void imposeB(typename Sys_eqs::stype num , grid_dist_iterator_sub<Sys_eqs::dims,typename Sys_eqs::b_grid::d_grid> it)
 	{
-#ifdef SE_CLASS1
-		if (Sys_eqs::num_cfields)
-			std::cerr << __FILE__ << ":" << __LINE__ << " if you do not provide ConstantFields in getMatrix, the system should not contain any constant field, while" << Sys_eqs::num_cfields << "\n";
-#endif
+		std::unordered_map<long int,float> cols;
+
+		// iterate all the grid points
+		while (it.isNext())
+		{
 
-		const ConstField fld[Sys_eqs::num_cfields];
+			b.get(row_b) = num;
 
-		getMatrix(g,fld);
+			++row_b;
+			++it;
+		}
 	}
 
 	/*! \brief produce the Matrix
@@ -90,31 +181,13 @@ public:
 	 *  \tparam Syst_eq System of equations, or a single equation
 	 *
 	 */
-	template<typename Grid> SparseMatrix<typename Sys_eqs::stype> getMatrix(const Grid & g, const ConstField(& fld)[Sys_eqs::num_cfields] )
+	template<typename Grid> SparseMatrix<typename Sys_eqs::stype> getMatrix()
 	{
-		// iterate across the full domain
-
-		auto it = g.getDomainIterator();
-		auto g_sm = g.getGrid();
-
-		Sys_eqs eqs();
-
-		// iterate all the grid points
-		while (it.isNext())
-		{
-			// get the position
-			auto key = it.get();
-
-			// convert into global coordinate the position
-			auto keyg = g.getGKey(key);
-
-			// Convert local key into global key
-			size_t row = g_sm.LidId(keyg);
+#ifdef SE_CLASS1
+		consistency();
+#endif
 
-			eqs.value(keyg,trpl);
 
-			++it;
-		}
 	}
 };
 

src/FiniteDifference/FDScheme_unit_tests.hpp

@@ -59,7 +59,7 @@ struct syss_nn
 	static const unsigned int nvar = 1;
 	static const unsigned int ord = EQS_FIELD;
 
-	static const unsigned int grid = STAGGERED_GRID;
+	static const unsigned int grid_type = STAGGERED_GRID;
 
 	// boundary at X and Y
 	static const bool boundary[];
@@ -78,7 +78,7 @@ struct syss_pp
 	static const unsigned int nvar = 1;
 	static const unsigned int ord = EQS_FIELD;
 
-	static const unsigned int grid = STAGGERED_GRID;
+	static const unsigned int grid_type = STAGGERED_GRID;
 
 	// boundary at X and Y
 	static const bool boundary[];
@@ -312,6 +312,52 @@ BOOST_AUTO_TEST_CASE( fd_test_use_periodic)
 	BOOST_REQUIRE_EQUAL(cols_y[14*16+15],-1);
 }
 
+BOOST_AUTO_TEST_CASE( fd_test_use_staggered_non_periodic)
+{
+	// grid size
+	size_t sz[2]={16,16};
+
+	// grid_sm
+	grid_sm<2,void> ginfo(sz);
+
+	// Create a derivative row Matrix
+	grid_key_dx<2> key11(1,1);
+	grid_key_dx<2> key00(0,0);
+	grid_key_dx<2> key22(2,2);
+	grid_key_dx<2> key1515(15,15);
+
+	// filled colums
+	std::unordered_map<long int,float> cols_x;
+	std::unordered_map<long int,float> cols_y;
+
+	D<x,Field<V,syss_pp>,syss_pp>::value(key11,ginfo,cols_x,1);
+	D<y,Field<V,syss_pp>,syss_pp>::value(key11,ginfo,cols_y,1);
+
+	BOOST_REQUIRE_EQUAL(cols_x.size(),2);
+	BOOST_REQUIRE_EQUAL(cols_y.size(),2);
+
+	BOOST_REQUIRE_EQUAL(cols_x[17+1],1);
+	BOOST_REQUIRE_EQUAL(cols_x[17],-1);
+
+	BOOST_REQUIRE_EQUAL(cols_y[17+16],1);
+	BOOST_REQUIRE_EQUAL(cols_y[17],-1);
+
+	cols_x.clear();
+	cols_y.clear();
+
+	D<x,Field<V,syss_pp>,syss_pp>::value(key00,ginfo,cols_x,1);
+	D<y,Field<V,syss_pp>,syss_pp>::value(key00,ginfo,cols_y,1);
+
+	BOOST_REQUIRE_EQUAL(cols_x.size(),2);
+	BOOST_REQUIRE_EQUAL(cols_y.size(),2);
+
+	BOOST_REQUIRE_EQUAL(cols_x[1],1);
+	BOOST_REQUIRE_EQUAL(cols_x[0],-1);
+
+	BOOST_REQUIRE_EQUAL(cols_y[16],1);
+	BOOST_REQUIRE_EQUAL(cols_y[0],-1);
+}
+
 /////////////// Laplacian test
 
 BOOST_AUTO_TEST_CASE( fd_test_lap_use_periodic)

src/FiniteDifference/eq.hpp

@@ -132,7 +132,7 @@ public:
 	static void value(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
 	{
 		if (Sys_eqs::ord == EQS_FIELD)
-			cols[gs.LinId(pos) + f] += coeff;
+			cols[gs.LinId(pos)*Sys_eqs::nvar + f] += coeff;
 		else
 			cols[gs.LinId(pos) + f * gs.size()] += coeff;
 	}

src/FiniteDifference/eq_unit_test.hpp

@@ -44,13 +44,16 @@ const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC};
 // Constant Field
 struct eta
 {
-	float val()	{return 1.0;}
+	typedef void const_field;
+
+	static float val()	{return 1.0;}
 };
 
 // Model the equations
 
 constexpr unsigned int v[] = {0,1};
 constexpr unsigned int P = 2;
+constexpr unsigned int ic = 2;
 
 typedef Field<v[x],lid_nn> v_x;
 typedef Field<v[y],lid_nn> v_y;
@@ -72,7 +75,8 @@ typedef sum<eta_lap_vy,p_x,lid_nn> vy_eq;
 
 typedef D<x,v_x,lid_nn> dx_vx;
 typedef D<y,v_y,lid_nn> dy_vy;
-typedef sum<dx_vx,dy_vy> incompressibility;
+typedef sum<dx_vx,dy_vy,lid_nn> ic_eq;
+
 
 BOOST_AUTO_TEST_SUITE( eq_test_suite )
 
@@ -87,6 +91,7 @@ BOOST_AUTO_TEST_CASE( lid_driven_cavity )
 	Ghost<2,float> g(0.01);
 
 	size_t sz[] = {32,32};
+	Padding<2> pd({1,2},{1,2});
 
 	// Initialize the global VCluster
 	init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
@@ -95,19 +100,79 @@ BOOST_AUTO_TEST_CASE( lid_driven_cavity )
 	grid_dist_id<2,float,scalar<float>,CartDecomposition<2,float>> g_dist(sz,domain,g);
 
 	// Distributed grid
-	FDScheme<lid_nn> fd;
+	FDScheme<lid_nn> fd(pd, g_dist.getGridInfo());
 
 	// start and end of the bulk
-	grid_key_dx<2> bulk_start(1,1);
-	grid_key_dx<2> bulk_end(sz[0]-1,sz[1]-1);
-
-	// Impose the vx equation
-	vx_eq vx;
-	fd.impose(vx, g_dist.getGridInfo() , g_dist.getSubDomainIterator(bulk_start,bulk_end));
-
-	// Impose the vy equation
-	vy_eq vy;
-	fd.impose(vy, g_dist.getGridInfo(), g_dist.getSubDomainIterator(bulk_start,bulk_end));
+	grid_key_dx<2> bulk_start(0,0);
+	grid_key_dx<2> bulk_end(sz[0],sz[1]);
+
+	// Impose the vx and vy equation in the bulk
+	fd.imposeA(vx_eq(), g_dist.getSubDomainIterator(bulk_start,bulk_end));
+	fd.imposeB(0.0, g_dist.getSubDomainIterator(bulk_start,bulk_end));
+	fd.imposeA(vy_eq(), g_dist.getSubDomainIterator(bulk_start,bulk_end));
+	fd.imposeB(0.0, g_dist.getSubDomainIterator(bulk_start,bulk_end));
+	fd.imposeA(ic_eq(), g_dist.getSubDomainIterator(bulk_start,bulk_end));
+	fd.imposeB(0.0, g_dist.getSubDomainIterator(bulk_start,bulk_end));
+
+	// Boundary condition on the left
+	grid_key_dx<2> bleft_start(0,-1);
+	grid_key_dx<2> bleft_end(0,33);
+	fd.imposeA(v_x(), g_dist.getSubDomainIterator(bleft_start,bleft_end));
+	fd.imposeB(0.0, g_dist.getSubDomainIterator(bleft_start,bleft_end));
+
+	// Boundary condition on the right
+	grid_key_dx<2> bright_start(32,-1);
+	grid_key_dx<2> bright_end(32,33);
+	fd.imposeA(v_x(), g_dist.getSubDomainIterator(bright_start,bright_end));
+	fd.imposeB(0.0, g_dist.getSubDomainIterator(bright_start,bright_end));
+
+	// Padding condition Pressure Left
+	grid_key_dx<2> pleft_start(-1,0);
+	grid_key_dx<2> pleft_end(-1,31);
+	fd.imposeA(Prs(), g_dist.getSubDomainIterator(bleft_start,bleft_end));
+	fd.imposeB(0.0, g_dist.getSubDomainIterator(bleft_start,bleft_end));
+
+	// Padding condition Pressure Right
+	grid_key_dx<2> pright_start(32,-1);
+	grid_key_dx<2> pright_end(33,33);
+	fd.imposeA(Prs(), g_dist.getSubDomainIterator(pright_start,pright_end));
+	fd.imposeB(0.0, g_dist.getSubDomainIterator(pright_start,pright_end));
+
+	// Padding condition vy Right
+	grid_key_dx<2> pvright_start(33,-1);
+	grid_key_dx<2> pvright_end(33,33);
+	fd.imposeA(v_y(), g_dist.getSubDomainIterator(pvright_start,pright_end));
+	fd.imposeB(0.0, g_dist.getSubDomainIterator(pvright_start,pright_end));
+
+	// Padding Bottom pressure
+	grid_key_dx<2> pbot_start(0,-1);
+	grid_key_dx<2> pbot_end(31,-1);
+	fd.imposeA(Prs(), g_dist.getSubDomainIterator(pbot_start,pbot_end));
+	fd.imposeB(0.0, g_dist.getSubDomainIterator(pbot_start,pbot_end));
+
+	// Bottom boundary condition vy
+	grid_key_dx<2> bbot_start(-1,0);
+	grid_key_dx<2> bbot_end(33,0);
+	fd.imposeA(v_y(), g_dist.getSubDomainIterator(bbot_start,bbot_end));
+	fd.imposeB(0.0, g_dist.getSubDomainIterator(bbot_start,bbot_end));
+
+	// Padding top pressure
+	grid_key_dx<2> ptop_start(0,32);
+	grid_key_dx<2> ptop_end(31,33);
+	fd.imposeA(Prs(), g_dist.getSubDomainIterator(ptop_start,ptop_end));
+	fd.imposeB(0.0, g_dist.getSubDomainIterator(ptop_start,ptop_end));
+
+	// Top boundary condition v_y
+	grid_key_dx<2> btop_start(-1,32);
+	grid_key_dx<2> btop_end(33,32);
+	fd.imposeA(v_y(), g_dist.getSubDomainIterator(btop_start,btop_end));
+	fd.imposeB(1.0, g_dist.getSubDomainIterator(btop_start,btop_end));
+
+	// Padding top v_x
+	grid_key_dx<2> pvtop_start(0,33);
+	grid_key_dx<2> pvtop_end(31,33);
+	fd.imposeA(v_x(), g_dist.getSubDomainIterator(pvtop_start,pvtop_end));
+	fd.imposeB(0.0, g_dist.getSubDomainIterator(pvtop_start,pvtop_end));
 }
 
 BOOST_AUTO_TEST_SUITE_END()

src/FiniteDifference/sum.hpp

@@ -36,7 +36,7 @@ struct sum_functor_value
 	/*! \brief constructor
 	 *
 	 */
-	sum_functor_value(grid_key_dx<last::dims> & key, const grid_sm<last::dims,void> & gs, typename last::stype coeff)
+	sum_functor_value(grid_key_dx<last::dims> & key, const grid_sm<last::dims,void> & gs, std::unordered_map<long int,typename last::stype> & cols, typename last::stype coeff)
 	:cols(cols),gs(gs),key(key),coeff(coeff)
 	{};
 
@@ -75,7 +75,7 @@ struct sum
 	inline static void value(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
 	{
 		// Sum functor
-		sum_functor_value<v_expr> sm(pos,gs,coeff);
+		sum_functor_value<v_expr> sm(pos,gs,cols,coeff);
 
 		// for each element in the expression calculate the non-zero Matrix elements
 		boost::mpl::for_each_ref< boost::mpl::range_c<int,0,v_sz::type::value - 1> >(sm);