OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
FD_simple_unit_test.cpp
1 //
2 // Created by jstark on 16.06.21.
3 //
4 
5 #define BOOST_TEST_DYN_LINK
6 #include <boost/test/unit_test.hpp>
7 
8 // Include header files for testing
9 #include "level_set/redistancing_Sussman/tests/l_norms/LNorms.hpp"
10 #include "Gaussian.hpp"
11 #include "FiniteDifference/FD_simple.hpp"
12 #include "level_set/redistancing_Sussman/HelpFunctions.hpp"
13 
14 BOOST_AUTO_TEST_SUITE(FDOrder1TestSuite)
15  const size_t Field = 0;
16  const size_t AnalyticalGradient = 1;
17  const size_t NumericalGradient = 2;
18  const size_t Error = 3;
19 
20 
21  const double EPSILON = std::numeric_limits<double>::epsilon();
22 
23  BOOST_AUTO_TEST_CASE(Forward_difference_1D_test)
24  {
25  double l2_norms [] = {0.573022, 0.288525, 0.171455};
26  const size_t grid_dim = 1;
27 
28  const double box_lower = -1.0;
29  const double box_upper = 1.0;
30  Box<grid_dim, double> box({box_lower}, {box_upper});
31  Ghost<grid_dim, long int> ghost(3);
32  typedef aggregate<double, Point<grid_dim, double>, Point<grid_dim, double>, double> props;
33  typedef grid_dist_id<grid_dim, double, props> grid_in_type;
34 
35  double mu = 0.5 * (box_upper - abs(box_lower));
36  double sigma = 0.1 * (box_upper - box_lower);
37 
38  int count = 0;
39  for (size_t N = 32; N <= 128; N *= 2, ++count)
40  {
41  const size_t sz[grid_dim] = {N};
42  grid_in_type g_dist(sz, box, ghost);
43 
44  g_dist.setPropNames({"Field", "AnalyticalGradient", "NumericalGradient", "Error"});
45 
46  auto gdom = g_dist.getDomainGhostIterator();
47  while (gdom.isNext())
48  {
49  auto key = gdom.get();
50  Point<grid_dim, double> p = g_dist.getPos(key);
51  // Initialize grid and ghost with gaussian function
52  g_dist.getProp<Field>(key) = gaussian(p, mu, sigma);
53  ++gdom;
54  }
55 
56  auto dom = g_dist.getDomainIterator();
57  while (dom.isNext())
58  {
59  auto key = dom.get();
60  Point<grid_dim, double> p = g_dist.getPos(key);
61 
62  for (int d = 0; d < grid_dim; d++)
63  {
64  // Analytical gradient
65  g_dist.getProp<AnalyticalGradient>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<Field>(key);
66  // 1st order Finite Difference gradient
67  g_dist.getProp<NumericalGradient>(key)[d] = FD_forward<Field>(g_dist, key, d);
68  }
69  ++dom;
70  }
71 
72  // Get the error between analytical and numerical solution
73 // get_absolute_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
74  get_relative_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
75 
76  LNorms<double> lNorms;
77  lNorms.get_l_norms_grid<Error>(g_dist);
78  BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms[count] + 0.00001 + EPSILON, "Checking L2-norm forward FD");
79 // write_lnorms_to_file(N, lNorms, "l_norms_FDfwd", "./");
80  std::cout << N << ", " << lNorms.l2 << ", " << lNorms.linf << std::endl;
81 // if (N==128) g_dist.write("grid_gaussian_FDfwd_N" + std::to_string(N), FORMAT_BINARY);
82  }
83  }
84  BOOST_AUTO_TEST_CASE(Backward_difference_1D_test)
85  {
86  double l2_norms [] = {0.573022, 0.288525, 0.171455};
87  const size_t grid_dim = 1;
88 
89  const double box_lower = -1.0;
90  const double box_upper = 1.0;
91  Box<grid_dim, double> box({box_lower}, {box_upper});
92  Ghost<grid_dim, long int> ghost(3);
93  typedef aggregate<double, Point<grid_dim, double>, Point<grid_dim, double>, double> props;
94  typedef grid_dist_id<grid_dim, double, props> grid_in_type;
95 
96  double mu = 0.5 * (box_upper - abs(box_lower));
97  double sigma = 0.1 * (box_upper - box_lower);
98 
99  int count = 0;
100  for (size_t N = 32; N <= 128; N *= 2, ++count)
101  {
102  const size_t sz[grid_dim] = {N};
103  grid_in_type g_dist(sz, box, ghost);
104 
105  g_dist.setPropNames({"Field", "AnalyticalGradient", "NumericalGradient", "Error"});
106 
107  auto gdom = g_dist.getDomainGhostIterator();
108  while (gdom.isNext())
109  {
110  auto key = gdom.get();
111  Point<grid_dim, double> p = g_dist.getPos(key);
112  // Initialize grid and ghost with gaussian function
113  g_dist.getProp<Field>(key) = gaussian(p, mu, sigma);
114  ++gdom;
115  }
116 
117  auto dom = g_dist.getDomainIterator();
118  while (dom.isNext())
119  {
120  auto key = dom.get();
121  Point<grid_dim, double> p = g_dist.getPos(key);
122 
123  for (int d = 0; d < grid_dim; d++)
124  {
125  // Analytical gradient
126  g_dist.getProp<AnalyticalGradient>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<Field>(key);
127  // 1st order Finite Difference gradient
128  g_dist.getProp<NumericalGradient>(key)[d] = FD_backward<Field>(g_dist, key, d);
129  }
130  ++dom;
131  }
132 
133  // Get the error between analytical and numerical solution
134 // get_absolute_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
135  get_relative_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
136 
137  LNorms<double> lNorms;
138  lNorms.get_l_norms_grid<Error>(g_dist);
139  BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms[count] + 0.00001 + EPSILON, "Checking L2-norm backward FD");
140 // write_lnorms_to_file(N, lNorms, "l_norms_FDbwd", "./");
141  std::cout << N << ", " << lNorms.l2 << ", " << lNorms.linf << std::endl;
142 // if (N==128) g_dist.write("grid_gaussian_FDbwd_N" + std::to_string(N), FORMAT_BINARY);
143  }
144  }
145  BOOST_AUTO_TEST_CASE(Central_difference_1D_test)
146  {
147  double l2_norms [] = {0.182302, 0.0405274, 0.00968203};
148  const size_t grid_dim = 1;
149 
150  const double box_lower = -1.0;
151  const double box_upper = 1.0;
152  Box<grid_dim, double> box({box_lower}, {box_upper});
153  Ghost<grid_dim, long int> ghost(3);
154  typedef aggregate<double, Point<grid_dim, double>, Point<grid_dim, double>, double> props;
155  typedef grid_dist_id<grid_dim, double, props> grid_in_type;
156 
157  double mu = 0.5 * (box_upper - abs(box_lower));
158  double sigma = 0.1 * (box_upper - box_lower);
159 
160  int count = 0;
161  for (size_t N = 32; N <= 128; N *= 2, ++count)
162  {
163  const size_t sz[grid_dim] = {N};
164  grid_in_type g_dist(sz, box, ghost);
165 
166  g_dist.setPropNames({"Field", "AnalyticalGradient", "NumericalGradient", "Error"});
167 
168  auto gdom = g_dist.getDomainGhostIterator();
169  while (gdom.isNext())
170  {
171  auto key = gdom.get();
172  Point<grid_dim, double> p = g_dist.getPos(key);
173  // Initialize grid and ghost with gaussian function
174  g_dist.getProp<Field>(key) = gaussian(p, mu, sigma);
175  ++gdom;
176  }
177 
178  auto dom = g_dist.getDomainIterator();
179  while (dom.isNext())
180  {
181  auto key = dom.get();
182  Point<grid_dim, double> p = g_dist.getPos(key);
183 
184  for (int d = 0; d < grid_dim; d++)
185  {
186  // Analytical gradient
187  g_dist.getProp<AnalyticalGradient>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<Field>(key);
188  // 1st order Finite Difference gradient
189  g_dist.getProp<NumericalGradient>(key)[d] = FD_central<Field>(g_dist, key, d);
190  }
191  ++dom;
192  }
193 
194  // Get the error between analytical and numerical solution
195 // get_absolute_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
196  get_relative_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
197 
198  LNorms<double> lNorms;
199  lNorms.get_l_norms_grid<Error>(g_dist);
200  BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms[count] + 0.00001 + EPSILON, "Checking L2-norm central FD");
201 // write_lnorms_to_file(N, lNorms, "l_norms_FDbwd", "./");
202  std::cout << N << ", " << lNorms.l2 << ", " << lNorms.linf << std::endl;
203 // if (N==128) g_dist.write("grid_gaussian_FDbwd_N" + std::to_string(N), FORMAT_BINARY);
204  }
205  }
206 BOOST_AUTO_TEST_SUITE_END()
Field
Definition: eq.hpp:82
LNorms.hpp
Header file containing functions for computing the error and the L_2 / L_infinity norm.
Point
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:27
Ghost
Definition: Ghost.hpp:39
grid_dist_id
This is a distributed grid.
Definition: grid_dist_id.hpp:96
Error
Out-of-bound policy kill the program.
Definition: vector_dist_ofb.hpp:69
Point::get
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition: Point.hpp:172
LNorms::get_l_norms_grid
void get_l_norms_grid(gridtype &grid)
Computes the L_2 and L_infinity norm on the basis of the precomputed error on a grid.
Definition: LNorms.hpp:123
Box
This class represent an N-dimensional box.
Definition: Box.hpp:60
aggregate
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition: aggregate.hpp:214
LNorms
Class for computing the l2/l_infinity norm for distributed grids and vectors based on given errors.
Definition: LNorms.hpp:103
HelpFunctions.hpp
Header file containing small mathematical help-functions that are needed e.g. for the Sussman redista...