OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
Upwind_gradient_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/Upwind_gradient.hpp"
12 #include "level_set/redistancing_Sussman/HelpFunctions.hpp"
13 
14 BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
15  const double EPSILON = 1e-14;
16  const size_t f_gaussian = 0;
17  const size_t Sign = 1;
18  const size_t df_gaussian = 2;
19  const size_t df_upwind = 3;
20  const size_t Error = 4;
21 
22  BOOST_AUTO_TEST_CASE(Upwind_gradient_order1_1D_test)
23  {
24  const int convergence_order = 1;
25  const size_t grid_dim = 1;
26 
27  const double box_lower = -1.0;
28  const double box_upper = 1.0;
29  Box<grid_dim, double> box({box_lower}, {box_upper});
30  Ghost<grid_dim, long int> ghost(3);
31  typedef aggregate<double, int, Point<grid_dim, double>, Point<grid_dim, double>, double> props;
32  typedef grid_dist_id<grid_dim, double, props> grid_in_type;
33 
34  double mu = 0.5 * (box_upper - abs(box_lower));
35  double sigma = 0.3 * (box_upper - box_lower);
36  int count = 0;
37  size_t N = 32;
38 // for (size_t N = 32; N <= 128; N *= 2, ++count)
39  {
40  const size_t sz[grid_dim] = {N};
41  grid_in_type g_dist(sz, box, ghost);
42 
43  g_dist.setPropNames({"f_gaussian", "Sign", "df_gaussian", "df_upwind", "Error"});
44 
45  auto gdom = g_dist.getDomainGhostIterator();
46  while (gdom.isNext())
47  {
48  auto key = gdom.get();
49  Point<grid_dim, double> p = g_dist.getPos(key);
50  // Initialize grid and ghost with gaussian function
51  g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
52  // Initialize the sign of the gaussian function
53  g_dist.getProp<Sign>(key) = sgn(g_dist.getProp<f_gaussian>(key));
54  ++gdom;
55  }
56 
57  auto dom = g_dist.getDomainIterator();
58  while (dom.isNext())
59  {
60  auto key = dom.get();
61  Point<grid_dim, double> p = g_dist.getPos(key);
62 
63  for (int d = 0; d < grid_dim; d++)
64  {
65  // Analytical gradient
66  g_dist.getProp<df_gaussian>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<f_gaussian>(key);
67  }
68  ++dom;
69  }
70 
71  // Upwind gradient with specific convergence order, not becoming one-sided at the boundary
72  get_upwind_gradient<f_gaussian, Sign, df_upwind>(g_dist, convergence_order, false);
73 
74  // Get the error between analytical and numerical solution
75  get_relative_error<df_upwind, df_gaussian, Error>(g_dist);
76 
77  LNorms<double> lNorms;
78  lNorms.get_l_norms_grid<Error>(g_dist);
79 
80  if (N==32) BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.38954184748195 + EPSILON, "Checking L2-norm upwind gradient "
81  "order " + std::to_string
82  (convergence_order));
83 
84 
85  }
86  }
87  BOOST_AUTO_TEST_CASE(Upwind_gradient_order3_1D_test)
88  {
89  const int convergence_order = 3;
90  const size_t grid_dim = 1;
91 
92  const double box_lower = -1.0;
93  const double box_upper = 1.0;
94  Box<grid_dim, double> box({box_lower}, {box_upper});
95  Ghost<grid_dim, long int> ghost(3);
96  typedef aggregate<double, int, Point<grid_dim, double>, Point<grid_dim, double>, double> props;
97  typedef grid_dist_id<grid_dim, double, props> grid_in_type;
98 
99  double mu = 0.5 * (box_upper - abs(box_lower));
100  double sigma = 0.3 * (box_upper - box_lower);
101  int count = 0;
102  size_t N = 32;
103 // for (size_t N = 32; N <= 128; N *= 2, ++count)
104  {
105  const size_t sz[grid_dim] = {N};
106  grid_in_type g_dist(sz, box, ghost);
107 
108  g_dist.setPropNames({"f_gaussian", "Sign", "df_gaussian", "df_upwind", "Error"});
109 
110  auto gdom = g_dist.getDomainGhostIterator();
111  while (gdom.isNext())
112  {
113  auto key = gdom.get();
114  Point<grid_dim, double> p = g_dist.getPos(key);
115  // Initialize grid and ghost with gaussian function
116  g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
117  // Initialize the sign of the gaussian function
118  g_dist.getProp<Sign>(key) = sgn(g_dist.getProp<f_gaussian>(key));
119  ++gdom;
120  }
121 
122  auto dom = g_dist.getDomainIterator();
123  while (dom.isNext())
124  {
125  auto key = dom.get();
126  Point<grid_dim, double> p = g_dist.getPos(key);
127 
128  for (int d = 0; d < grid_dim; d++)
129  {
130  // Analytical gradient
131  g_dist.getProp<df_gaussian>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<f_gaussian>(key);
132  }
133  ++dom;
134  }
135 
136  // Upwind gradient with specific convergence order, not becoming one-sided at the boundary
137  get_upwind_gradient<f_gaussian, Sign, df_upwind>(g_dist, convergence_order, false);
138 
139  // Get the error between analytical and numerical solution
140  get_relative_error<df_upwind, df_gaussian, Error>(g_dist);
141 
142  LNorms<double> lNorms;
143  lNorms.get_l_norms_grid<Error>(g_dist);
144 
145  if (N==32) BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.08667855716144 + EPSILON, "Checking L2-norm upwind gradient "
146  "order "
147  + std::to_string(convergence_order));
148  }
149  }
150  BOOST_AUTO_TEST_CASE(Upwind_gradient_order5_1D_test)
151  {
152  const int convergence_order = 5;
153  const size_t grid_dim = 1;
154 
155  const double box_lower = -1.0;
156  const double box_upper = 1.0;
157  Box<grid_dim, double> box({box_lower}, {box_upper});
158  Ghost<grid_dim, long int> ghost(3);
159  typedef aggregate<double, int, Point<grid_dim, double>, Point<grid_dim, double>, double> props;
160  typedef grid_dist_id<grid_dim, double, props> grid_in_type;
161 
162  double mu = 0.5 * (box_upper - abs(box_lower));
163  double sigma = 0.3 * (box_upper - box_lower);
164  int count = 0;
165  size_t N = 32;
166 // for (size_t N = 32; N <= 128; N *= 2, ++count)
167  {
168  const size_t sz[grid_dim] = {N};
169  grid_in_type g_dist(sz, box, ghost);
170 
171  g_dist.setPropNames({"f_gaussian", "Sign", "df_gaussian", "df_upwind", "Error"});
172 
173  auto gdom = g_dist.getDomainGhostIterator();
174  while (gdom.isNext())
175  {
176  auto key = gdom.get();
177  Point<grid_dim, double> p = g_dist.getPos(key);
178  // Initialize grid and ghost with gaussian function
179  g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
180  // Initialize the sign of the gaussian function
181  g_dist.getProp<Sign>(key) = sgn(g_dist.getProp<f_gaussian>(key));
182  ++gdom;
183  }
184 
185  auto dom = g_dist.getDomainIterator();
186  while (dom.isNext())
187  {
188  auto key = dom.get();
189  Point<grid_dim, double> p = g_dist.getPos(key);
190 
191  for (int d = 0; d < grid_dim; d++)
192  {
193  // Analytical gradient
194  g_dist.getProp<df_gaussian>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<f_gaussian>(key);
195  }
196  ++dom;
197  }
198 
199  // Upwind gradient with specific convergence order, not becoming one-sided at the boundary
200  get_upwind_gradient<f_gaussian, Sign, df_upwind>(g_dist, convergence_order, false);
201 
202  // Get the error between analytical and numerical solution
203  get_relative_error<df_upwind, df_gaussian, Error>(g_dist);
204 
205  LNorms<double> lNorms;
206  lNorms.get_l_norms_grid<Error>(g_dist);
207 
208  if (N==32) BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.03215172234342 + EPSILON, "Checking L2-norm upwind gradient "
209  "order "
210  + std::to_string(convergence_order));
211  }
212  }
213  BOOST_AUTO_TEST_CASE(Upwind_gradient_3D_test)
214  {
215  {
216  const size_t grid_dim = 3;
217  const double box_lower = -1.0;
218  const double box_upper = 1.0;
219  Box<grid_dim, double> box({box_lower, box_lower, box_lower}, {box_upper, box_upper, box_upper});
220  Ghost<grid_dim, long int> ghost(3);
221  typedef aggregate<double, int, Point<grid_dim, double>, Point<grid_dim, double>, double> props;
222  typedef grid_dist_id<grid_dim, double, props> grid_in_type;
223 
224  double mu = 0.5 * (box_upper - abs(box_lower));
225  double sigma = 0.3 * (box_upper - box_lower);
226 
227  size_t N = 32;
228 // for(size_t N = 32; N <=256; N *=2)
229  {
230  const size_t sz[grid_dim] = {N, N, N};
231  grid_in_type g_dist(sz, box, ghost);
232  g_dist.setPropNames({"f_gaussian", "Sign", "df_gaussian", "df_upwind", "Error"});
233 
234  auto gdom = g_dist.getDomainGhostIterator();
235  while (gdom.isNext())
236  {
237  auto key = gdom.get();
238  Point<grid_dim, double> p = g_dist.getPos(key);
239  // Initialize grid and ghost with gaussian function
240  g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
241  // Initialize the sign of the gaussian function
242  g_dist.getProp<Sign>(key) = sgn(g_dist.getProp<f_gaussian>(key));
243  ++gdom;
244  }
245 
246  auto dom = g_dist.getDomainIterator();
247  while (dom.isNext())
248  {
249  auto key = dom.get();
250  Point<grid_dim, double> p = g_dist.getPos(key);
251 
252  for (int d = 0; d < grid_dim; d++)
253  {
254  // Analytical gradient
255  g_dist.getProp<df_gaussian>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<f_gaussian>(key);
256  }
257  ++dom;
258  }
259 
260  for (int convergence_order = 1; convergence_order <=5; convergence_order +=2 )
261  {
262  // Upwind gradient with specific convergence order, not becoming one-sided at the boundary
263  get_upwind_gradient<f_gaussian, Sign, df_upwind>(g_dist, convergence_order, false);
264  // Get the error between analytical and numerical solution
265  get_relative_error<df_upwind, df_gaussian, Error>(g_dist);
266 
267  LNorms<double> lNorms;
268  lNorms.get_l_norms_grid<Error>(g_dist);
269 
270  if(N==32)
271  {
272  switch(convergence_order)
273  {
274  case 1:
275  BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.38954184748194 + EPSILON, "Checking L2-norm upwind gradient order " + std::to_string(convergence_order));
276  break;
277  case 3:
278  BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.08667855716144 + EPSILON, "Checking L2-norm upwind gradient order " + std::to_string(convergence_order));
279  break;
280  case 5:
281  BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.02285996528578 + EPSILON, "Checking L2-norm upwind gradient order " + std::to_string(convergence_order));
282 
283  break;
284  default:
285  std::cout << "Checking only implemented for convergence order 1, 3 and 5." << std::endl;
286  }
287  }
288  }
289 
290  }
291  }
292  }
293 BOOST_AUTO_TEST_SUITE_END()
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
sgn
int sgn(T val)
Gets the sign of a variable.
Definition: HelpFunctions.hpp:24
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
Upwind_gradient.hpp
Approximating upwind gradients on a grid with the following options for the order of accuracy: 1,...
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...