OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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
18BOOST_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 }
316BOOST_AUTO_TEST_SUITE_END()
LNorms.hpp
Header file containing functions for computing the error and the L_2 / L_infinity norm.
NarrowBand.hpp
Class for getting the narrow band around the interface.
PathsAndFiles.hpp
Header file containing functions for creating files and folders.
RedistancingSussman.hpp
Class for reinitializing a level-set function into a signed distance function using Sussman redistanc...
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
LNorms
Class for computing the l2/l_infinity norm for distributed grids and vectors based on given errors.
Definition LNorms.hpp:105
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:124
Point
This class implement the point shape in an N-dimensional space.
Definition Point.hpp:28
Point::get
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition Point.hpp:172
grid_dist_id
This is a distributed grid.
Definition grid_dist_id.hpp:278
aggregate
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition aggregate.hpp:215