OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
Eno_Weno_unit_test.cpp
1 //
2 // Created by jstark on 14.06.21.
3 //
4 
5 #define BOOST_TEST_DYN_LINK
6 //#define BOOST_TEST_MAIN // in only one cpp file
7 #include <boost/test/unit_test.hpp>
8 
9 // Include redistancing files
10 #include "util/PathsAndFiles.hpp"
11 #include "level_set/redistancing_Sussman/RedistancingSussman.hpp"
12 #include "level_set/redistancing_Sussman/NarrowBand.hpp"
13 // Include header files for testing
14 #include "Draw/DrawSphere.hpp"
15 #include "level_set/redistancing_Sussman/tests/l_norms/LNorms.hpp"
16 #include "Gaussian.hpp"
17 
18 BOOST_AUTO_TEST_SUITE(EnoWenoTestSuite)
19  const size_t f_gaussian = 0;
20  const size_t df_gaussian = 1;
21  const size_t ENO_plus = 2;
22  const size_t ENO_minus = 3;
23  const size_t Error_plus = 4;
24  const size_t Error_minus = 5;
25 
26  double l2_norms_ENO [] = {0.084908, 0.014937, 0.002310};
27  double linf_norms_ENO [] = {0.228915, 0.047215, 0.008046};
28 
29  double l2_norms_WENO [] = {0.032770, 0.000935, 0.000052};
30  double linf_norms_WENO [] = {0.168860, 0.003369, 0.000232};
31 
32  BOOST_AUTO_TEST_CASE(Eno_1D_1stDerivative_test)
33  {
34  int count = 0;
35  for (size_t N = 32; N <= 128; N *= 2)
36  {
37  const size_t grid_dim = 1;
38  const size_t sz[grid_dim] = {N};
39  const double box_lower = -1.0;
40  const double box_upper = 1.0;
41  Box<grid_dim, double> box({box_lower}, {box_upper});
42  Ghost<grid_dim, long int> ghost(3);
43  typedef aggregate<double, Point<grid_dim, double>, Point<grid_dim, double>, Point<grid_dim, double>, double, double> props;
44  typedef grid_dist_id<grid_dim, double, props> grid_in_type;
45  grid_in_type g_dist(sz, box, ghost);
46 
47  g_dist.setPropNames({"f_gaussian", "df_gaussian", "ENO_plus", "ENO_minus", "Error_plus", "Error_minus"});
48 
49  double mu = 0.5 * (box_upper - abs(box_lower));
50  double sigma = 0.3 * (box_upper - box_lower);
51 
52  auto gdom = g_dist.getDomainGhostIterator();
53  while (gdom.isNext())
54  {
55  auto key = gdom.get();
56  Point<grid_dim, double> p = g_dist.getPos(key);
57  // Initialize grid and ghost with gaussian function
58  g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
59  ++gdom;
60  }
61 
62  auto dom = g_dist.getDomainIterator();
63  while (dom.isNext())
64  {
65  auto key = dom.get();
66  Point<grid_dim, double> p = g_dist.getPos(key);
67 
68  for (int d = 0; d < grid_dim; d++)
69  {
70  // Analytical solution
71  g_dist.getProp<df_gaussian>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<f_gaussian>(key);
72  // ENO_plus
73  g_dist.getProp<ENO_plus>(key)[d] = ENO_3_Plus<f_gaussian>(g_dist, key, d);
74  // ENO_minus
75  g_dist.getProp<ENO_minus>(key)[d] = ENO_3_Minus<f_gaussian>(g_dist, key, d);
76  }
77 
78  ++dom;
79  }
80  // Get the error between analytical and numerical solution
81  get_relative_error<df_gaussian, ENO_plus, Error_plus>(g_dist);
82  get_relative_error<df_gaussian, ENO_minus, Error_minus>(g_dist);
83 
84  {
85  LNorms<double> lNorms;
86  lNorms.get_l_norms_grid<Error_plus>(g_dist);
87  BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms_ENO[count] + 0.000001, "Checking L2-norm ENO");
88  BOOST_CHECK_MESSAGE(lNorms.linf < linf_norms_ENO[count] + 0.000001, "Checking Linf-norm "
89  "ENO");
90 // write_lnorms_to_file(N, lNorms, "l_norms_ENO_plus", "./");
91  }
92 
93  {
94  LNorms<double> lNorms;
95  lNorms.get_l_norms_grid<Error_minus>(g_dist);
96  BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms_ENO[count] + 0.000001, "Checking L2-norm ENO");
97  BOOST_CHECK_MESSAGE(lNorms.linf < linf_norms_ENO[count] + 0.000001, "Checking Linf-norm "
98  "ENO");
99  }
100 // g_dist.write("grid_gaussian_ENO_1D_N" + std::to_string(N), FORMAT_BINARY);
101  count++;
102  }
103  }
104  const size_t WENO_plus = 2;
105  const size_t WENO_minus = 3;
106  BOOST_AUTO_TEST_CASE(Weno_1D_1stDerivative_test)
107  {
108  int count = 0;
109  for (size_t N = 32; N <= 128; N *= 2)
110  {
111  const size_t grid_dim = 1;
112  const size_t sz[grid_dim] = {N};
113  const double box_lower = -1.0;
114  const double box_upper = 1.0;
115  Box<grid_dim, double> box({box_lower}, {box_upper});
116  Ghost<grid_dim, long int> ghost(3);
117  typedef aggregate<double, Point<grid_dim, double>, Point<grid_dim, double>, Point<grid_dim, double>, double, double> props;
118  typedef grid_dist_id<grid_dim, double, props> grid_in_type;
119  grid_in_type g_dist(sz, box, ghost);
120 
121  g_dist.setPropNames({"f_gaussian", "df_gaussian", "WENO_plus", "WENO_minus", "Error_plus", "Error_minus"});
122 
123  double mu = 0.5 * (box_upper - abs(box_lower));
124  double sigma = 0.3 * (box_upper - box_lower);
125 
126  auto gdom = g_dist.getDomainGhostIterator();
127  while (gdom.isNext())
128  {
129  auto key = gdom.get();
130  Point<grid_dim, double> p = g_dist.getPos(key);
131  // Initialize grid and ghost with gaussian function
132  g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
133  ++gdom;
134  }
135 
136  auto dom = g_dist.getDomainIterator();
137  while (dom.isNext())
138  {
139  auto key = dom.get();
140  Point<grid_dim, double> p = g_dist.getPos(key);
141  for (int d = 0; d < grid_dim; d++)
142  {
143  // Analytical solution
144  g_dist.getProp<df_gaussian>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<f_gaussian>(key);
145  // WENO_plus
146  g_dist.getProp<WENO_plus>(key)[d] = WENO_5_Plus<f_gaussian>(g_dist, key, d);
147  // WENO_minus
148  g_dist.getProp<WENO_minus>(key)[d] = WENO_5_Minus<f_gaussian>(g_dist, key, d);
149  }
150 
151  ++dom;
152  }
153  // Get the error between analytical and numerical solution
154  get_relative_error<df_gaussian, WENO_plus, Error_plus>(g_dist);
155  get_relative_error<df_gaussian, WENO_minus, Error_minus>(g_dist);
156 
157  {
158  LNorms<double> lNorms;
159  lNorms.get_l_norms_grid<Error_plus>(g_dist);
160  BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms_WENO[count] + 0.000001, "Checking L2-norm WENO");
161  BOOST_CHECK_MESSAGE(lNorms.linf < linf_norms_WENO[count] + 0.000001, "Checking Linf-norm "
162  "WENO");
163 // write_lnorms_to_file(N, lNorms, "l_norms_WENO_plus", "./");
164  }
165 
166  {
167  LNorms<double> lNorms;
168  lNorms.get_l_norms_grid<Error_minus>(g_dist);
169  BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms_WENO[count] + 0.000001, "Checking L2-norm WENO");
170  BOOST_CHECK_MESSAGE(lNorms.linf < linf_norms_WENO[count] + 0.000001, "Checking Linf-norm "
171  "WENO");
172 // write_lnorms_to_file(N, lNorms, "l_norms_WENO_minus", "./");
173  }
174 // g_dist.write("grid_gaussian_WENO_1D_N" + std::to_string(N), FORMAT_BINARY);
175  count++;
176  }
177  }
178 
179  BOOST_AUTO_TEST_CASE(Eno_3D_test)
180  {
181  int count = 0;
182  for(size_t N=32; N<=128; N*=2)
183  {
184  const size_t grid_dim = 3;
185  const size_t sz[grid_dim] = {N, N, N};
186  const double box_lower = -1.0;
187  const double box_upper = 1.0;
188  Box<grid_dim, double> box({box_lower, box_lower, box_lower}, {box_upper, box_upper, box_upper});
189  Ghost<grid_dim, long int> ghost(3);
190  typedef aggregate<double, Point<grid_dim, double>, Point<grid_dim, double>, Point<grid_dim, double>, double, double> props;
191  typedef grid_dist_id<grid_dim, double, props> grid_in_type;
192  grid_in_type g_dist(sz, box, ghost);
193 
194  g_dist.setPropNames({"f_gaussian", "df_gaussian", "ENO_plus", "ENO_minus", "Error_plus", "Error_minus"});
195 
196  double mu = 0.5 * (box_upper - abs(box_lower));
197  double sigma = 0.3 * (box_upper - box_lower);
198 
199  auto gdom = g_dist.getDomainGhostIterator();
200  while (gdom.isNext())
201  {
202  auto key = gdom.get();
203  Point<grid_dim, double> p = g_dist.getPos(key);
204  // Initialize grid and ghost with gaussian function
205  g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
206  ++gdom;
207  }
208 
209  auto dom = g_dist.getDomainIterator();
210  while (dom.isNext())
211  {
212  auto key = dom.get();
213  Point<grid_dim, double> p = g_dist.getPos(key);
214  for (int d = 0; d < grid_dim; d++)
215  {
216  g_dist.getProp<df_gaussian>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<f_gaussian>(key);
217  g_dist.getProp<ENO_plus>(key)[d] = ENO_3_Plus<f_gaussian>(g_dist, key, d);
218  g_dist.getProp<ENO_minus>(key)[d] = ENO_3_Minus<f_gaussian>(g_dist, key, d);
219  }
220 
221  ++dom;
222  }
223 
224  // Get the error between analytical and numerical solution
225  get_relative_error<df_gaussian, ENO_plus, Error_plus>(g_dist);
226  get_relative_error<df_gaussian, ENO_minus, Error_minus>(g_dist);
227 
228  {
229  LNorms<double> lNorms;
230  lNorms.get_l_norms_grid<Error_plus>(g_dist);
231  BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms_ENO[count] + 0.000001, "Checking L2-norm ENO");
232  BOOST_CHECK_MESSAGE(lNorms.linf < linf_norms_ENO[count] + 0.000001, "Checking Linf-norm "
233  "ENO");
234 // write_lnorms_to_file(N, lNorms, "l_norms_3D_ENO_plus", "./");
235  }
236 
237  {
238  LNorms<double> lNorms;
239  lNorms.get_l_norms_grid<Error_minus>(g_dist);
240  BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms_ENO[count] + 0.000001, "Checking L2-norm ENO");
241  BOOST_CHECK_MESSAGE(lNorms.linf < linf_norms_ENO[count] + 0.000001, "Checking Linf-norm "
242  "ENO");
243  }
244 // g_dist.write("grid_gaussian_3D_N" + std::to_string(N), FORMAT_BINARY);
245  count++;
246  }
247  }
248  BOOST_AUTO_TEST_CASE(Weno_3D_test)
249  {
250  int count = 0;
251  for(size_t N=32; N<=128; N*=2)
252  {
253  const size_t grid_dim = 3;
254  const size_t sz[grid_dim] = {N, N, N};
255  const double box_lower = -1.0;
256  const double box_upper = 1.0;
257  Box<grid_dim, double> box({box_lower, box_lower, box_lower}, {box_upper, box_upper, box_upper});
258  Ghost<grid_dim, long int> ghost(3);
259  typedef aggregate<double, Point<grid_dim, double>, Point<grid_dim, double>, Point<grid_dim, double>, double, double> props;
260  typedef grid_dist_id<grid_dim, double, props> grid_in_type;
261  grid_in_type g_dist(sz, box, ghost);
262 
263  g_dist.setPropNames({"f_gaussian", "df_gaussian", "WENO_plus", "WENO_minus", "Error_plus", "Error_minus"});
264 
265  double mu = 0.5 * (box_upper - abs(box_lower));
266  double sigma = 0.3 * (box_upper - box_lower);
267 
268  auto gdom = g_dist.getDomainGhostIterator();
269  while (gdom.isNext())
270  {
271  auto key = gdom.get();
272  Point<grid_dim, double> p = g_dist.getPos(key);
273  // Initialize grid and ghost with gaussian function
274  g_dist.getProp<f_gaussian>(key) = gaussian(p, mu, sigma);
275  ++gdom;
276  }
277 
278  auto dom = g_dist.getDomainIterator();
279  while (dom.isNext())
280  {
281  auto key = dom.get();
282  Point<grid_dim, double> p = g_dist.getPos(key);
283  for (int d = 0; d < grid_dim; d++)
284  {
285  g_dist.getProp<df_gaussian>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<f_gaussian>(key);
286  g_dist.getProp<WENO_plus>(key)[d] = WENO_5_Plus<f_gaussian>(g_dist, key, d);
287  g_dist.getProp<WENO_minus>(key)[d] = WENO_5_Minus<f_gaussian>(g_dist, key, d);
288  }
289 
290  ++dom;
291  }
292 
293  // Get the error between analytical and numerical solution
294  get_relative_error<df_gaussian, WENO_plus, Error_plus>(g_dist);
295  get_relative_error<df_gaussian, WENO_minus, Error_minus>(g_dist);
296 
297  {
298  LNorms<double> lNorms;
299  lNorms.get_l_norms_grid<Error_plus>(g_dist);
300  BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms_WENO[count] + 0.000001, "Checking L2-norm WENO");
301  BOOST_CHECK_MESSAGE(lNorms.linf < linf_norms_WENO[count] + 0.000001, "Checking Linf-norm "
302  "WENO");
303  }
304 
305  {
306  LNorms<double> lNorms;
307  lNorms.get_l_norms_grid<Error_minus>(g_dist);
308  BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms_WENO[count] + 0.000001, "Checking L2-norm WENO");
309  BOOST_CHECK_MESSAGE(lNorms.linf < linf_norms_WENO[count] + 0.000001, "Checking Linf-norm "
310  "WENO");
311  }
312 // g_dist.write("grid_gaussian_3D_N" + std::to_string(N), FORMAT_BINARY);
313  count++;
314  }
315  }
316 BOOST_AUTO_TEST_SUITE_END()
RedistancingSussman.hpp
Class for reinitializing a level-set function into a signed distance function using Sussman redistanc...
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
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
PathsAndFiles.hpp
Header file containing functions for creating files and folders.
aggregate
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition: aggregate.hpp:214
NarrowBand.hpp
Class for getting the narrow band around the interface.
LNorms
Class for computing the l2/l_infinity norm for distributed grids and vectors based on given errors.
Definition: LNorms.hpp:103