Added central finite difference incl unit test.

src/FiniteDifference/FD_order1.hpp → src/FiniteDifference/FD_simple.hpp

 //
 // Created by jstark on 17.05.21.
 //
-
-#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FD_ORDER1_HPP
-#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FD_ORDER1_HPP
+#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FD_SIMPLE_HPP
+#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FD_SIMPLE_HPP
 
 // Include OpenFPM header files
 #include "Grid/grid_dist_id.hpp"
@@ -40,10 +39,78 @@ double FD_backward(gridtype & grid, keytype & key, size_t d)
 }
 
 
+/**@brief Computes the central finite difference of a scalar field on the current grid node.
+ *
+ * @tparam Field Size_t index of property for which the gradient should be computed.
+ * @tparam gridtype Type of input grid.
+ * @tparam keytype Type of key variable.
+ * @param grid Grid, on which the gradient should be computed.
+ * @param key Key that contains the index of the current grid node.
+ * @param d Variable (size_t) that contains the dimension.
+ * @return Partial central finite difference for dimension d of Field on the current node with index key.
+ */
+template <size_t Field, typename gridtype, typename keytype>
+double FD_central(gridtype & grid, keytype & key, size_t d)
+{
+	return (grid.template get<Field>(key.move(d, 1))
+			- grid.template get<Field>(key.move(d, -1)))
+			/ (2 * grid.getSpacing()[d]);
+}
 
+/**@brief Computes the central finite difference of a scalar field on the full grid.
+ * 
+ * @tparam Field Size_t index of input property for which the gradient should be computed (scalar field).
+ * @tparam Gradient Size_t index of output property (vector field).
+ * @tparam gridtype Type of input grid.
+ * @param grid Grid, on which the gradient should be computed.
+ * @param one_sided_BC Boolian. If true, stencil becoming 1st order one sided at the grid boundary (fwd/bwd FD). If 
+ * false, stencil remains centered 2nd order at the boundary, using the ghost layer.
+ */
+template <size_t Field, size_t Gradient, typename gridtype>
+void get_central_FD_grid(gridtype & grid, const bool one_sided_BC)
+{
+	grid.template ghost_get<Field>(KEEP_PROPERTIES);
+	auto dom = grid.getDomainIterator();
+	if (one_sided_BC)
+	{
+		while (dom.isNext())
+		{
+			auto key = dom.get();
+			auto key_g = grid.getGKey(key);
+			for(size_t d = 0; d < gridtype::dims; d++)
+			{
+				// Grid nodes inside and distance from boundary > stencil-width
+				if (key_g.get(d) > 0 && key_g.get(d) < grid.size(d) - 1) // if point lays with min. 1 nodes distance to
+					// boundary
+				{
+					grid.template get<Gradient> (key) [d] = FD_central<Field>(grid, key, d);
+				}
+				else if (key_g.get(d) == 0) // if point lays at left boundary, use right sided kernel
+				{
+					grid.template get<Gradient> (key) [d] = FD_forward<Field>(grid, key, d);
+				}
+				else if (key_g.get(d) >= grid.size(d) - 1) // if point lays at right boundary, use left sided kernel
+				{
+					grid.template get<Gradient> (key) [d] = FD_backward<Field>(grid, key, d);
+				}
+			}
+			++dom;
+		}
+	}
+	else
+	{
+		while (dom.isNext())
+		{
+			auto key = dom.get();
+			for(size_t d = 0; d < gridtype::dims; d++)
+			{
+				grid.template get<Gradient> (key) [d] = FD_central<Field>(grid, key, d);
+			}
+			++dom;
+		}
+	}
+}
 
 
-
-
-#endif //OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FD_ORDER1_HPP
+#endif //OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FD_SIMPLE_HPP
 

src/FiniteDifference/Upwind_gradient.hpp

@@ -25,7 +25,7 @@
 
 // Include OpenFPM header files
 #include "Grid/grid_dist_id.hpp"
-#include "FD_order1.hpp"
+#include "FD_simple.hpp"
 #include "Eno_Weno.hpp"
 
 /**@brief Upwinding: For a specific dimension, from the forward and backward gradient find the upwind side.

src/FiniteDifference/tests/FD_order1_unit_test.cpp → src/FiniteDifference/tests/FD_simple_unit_test.cpp

@@ -8,7 +8,7 @@
 // Include header files for testing
 #include "level_set/redistancing_Sussman/tests/l_norms/LNorms.hpp"
 #include "Gaussian.hpp"
-#include "FiniteDifference/Upwind_gradient.hpp"
+#include "FiniteDifference/FD_simple.hpp"
 #include "level_set/redistancing_Sussman/HelpFunctions.hpp"
 
 BOOST_AUTO_TEST_SUITE(FDOrder1TestSuite)
@@ -17,15 +17,12 @@ BOOST_AUTO_TEST_SUITE(FDOrder1TestSuite)
 	const size_t NumericalGradient  = 2;
 	const size_t Error              = 3;
 	
-	// 32, 0.573022, 1.33305
-	// 64, 0.288525, 1
-	//128, 0.171455, 1
-	double l2_norms [] = {0.573022, 0.288525, 0.171455};
 	
 	const double EPSILON = std::numeric_limits<double>::epsilon();
 	
 	BOOST_AUTO_TEST_CASE(Forward_difference_1D_test)
 	{
+		double l2_norms [] = {0.573022, 0.288525, 0.171455};
 		const size_t grid_dim  = 1;
 		
 		const double box_lower = -1.0;
@@ -78,7 +75,7 @@ BOOST_AUTO_TEST_SUITE(FDOrder1TestSuite)
 			
 			L_norms lNorms;
 			lNorms = get_l_norms_grid<Error>(g_dist);
-			BOOST_CHECK_MESSAGE(lNorms.l2   < l2_norms[count] + 0.00001 + EPSILON, "Checking L2-norm ENO");
+			BOOST_CHECK_MESSAGE(lNorms.l2   < l2_norms[count] + 0.00001 + EPSILON, "Checking L2-norm forward FD");
 //			write_lnorms_to_file(N, lNorms, "l_norms_FDfwd", "./");
 			std::cout << N << ", " << lNorms.l2 << ", " << lNorms.linf << std::endl;
 //			if (N==128) g_dist.write("grid_gaussian_FDfwd_N" + std::to_string(N), FORMAT_BINARY);
@@ -86,6 +83,7 @@ BOOST_AUTO_TEST_SUITE(FDOrder1TestSuite)
 	}
 	BOOST_AUTO_TEST_CASE(Backward_difference_1D_test)
 	{
+		double l2_norms [] = {0.573022, 0.288525, 0.171455};
 		const size_t grid_dim  = 1;
 		
 		const double box_lower = -1.0;
@@ -138,7 +136,68 @@ BOOST_AUTO_TEST_SUITE(FDOrder1TestSuite)
 			
 			L_norms lNorms;
 			lNorms = get_l_norms_grid<Error>(g_dist);
-			BOOST_CHECK_MESSAGE(lNorms.l2   < l2_norms[count] + 0.00001 + EPSILON, "Checking L2-norm ENO");
+			BOOST_CHECK_MESSAGE(lNorms.l2   < l2_norms[count] + 0.00001 + EPSILON, "Checking L2-norm backward FD");
+//			write_lnorms_to_file(N, lNorms, "l_norms_FDbwd", "./");
+			std::cout << N << ", " << lNorms.l2 << ", " << lNorms.linf << std::endl;
+//			if (N==128) g_dist.write("grid_gaussian_FDbwd_N" + std::to_string(N), FORMAT_BINARY);
+		}
+	}
+	BOOST_AUTO_TEST_CASE(Central_difference_1D_test)
+	{
+		double l2_norms [] = {0.182302, 0.0405274, 0.00968203};
+		const size_t grid_dim  = 1;
+		
+		const double box_lower = -1.0;
+		const double box_upper = 1.0;
+		Box<grid_dim, double> box({box_lower}, {box_upper});
+		Ghost<grid_dim, long int> ghost(3);
+		typedef aggregate<double, Point<grid_dim, double>, Point<grid_dim, double>, double> props;
+		typedef grid_dist_id<grid_dim, double, props> grid_in_type;
+		
+		double mu = 0.5 * (box_upper - abs(box_lower));
+		double sigma = 0.1 * (box_upper - box_lower);
+		
+		int count = 0;
+		for (size_t N = 32; N <= 128; N *= 2, ++count)
+		{
+			const size_t sz[grid_dim] = {N};
+			grid_in_type g_dist(sz, box, ghost);
+			
+			g_dist.setPropNames({"Field", "AnalyticalGradient", "NumericalGradient", "Error"});
+			
+			auto gdom = g_dist.getDomainGhostIterator();
+			while (gdom.isNext())
+			{
+				auto key = gdom.get();
+				Point<grid_dim, double> p = g_dist.getPos(key);
+				// Initialize grid and ghost with gaussian function
+				g_dist.getProp<Field>(key) = gaussian(p, mu, sigma);
+				++gdom;
+			}
+			
+			auto dom = g_dist.getDomainIterator();
+			while (dom.isNext())
+			{
+				auto key = dom.get();
+				Point<grid_dim, double> p = g_dist.getPos(key);
+				
+				for (int d = 0; d < grid_dim; d++)
+				{
+					// Analytical gradient
+					g_dist.getProp<AnalyticalGradient>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<Field>(key);
+					// 1st order Finite Difference gradient
+					g_dist.getProp<NumericalGradient>(key)[d]  = FD_central<Field>(g_dist, key, d);
+				}
+				++dom;
+			}
+			
+			// Get the error between analytical and numerical solution
+//			get_absolute_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
+			get_relative_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
+			
+			L_norms lNorms;
+			lNorms = get_l_norms_grid<Error>(g_dist);
+			BOOST_CHECK_MESSAGE(lNorms.l2   < l2_norms[count] + 0.00001 + EPSILON, "Checking L2-norm central FD");
 //			write_lnorms_to_file(N, lNorms, "l_norms_FDbwd", "./");
 			std::cout << N << ", " << lNorms.l2 << ", " << lNorms.linf << std::endl;
 //			if (N==128) g_dist.write("grid_gaussian_FDbwd_N" + std::to_string(N), FORMAT_BINARY);