OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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
10#include "Gaussian.hpp"
13
14BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
15 const double EPSILON = 1e-14;
16 const size_t F_GAUSSIAN = 0;
17 const size_t VELOCITY = 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});
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", "VELOCITY", "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<VELOCITY>(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, VELOCITY, 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});
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", "VELOCITY", "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<VELOCITY>(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, VELOCITY, 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});
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", "VELOCITY", "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<VELOCITY>(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, VELOCITY, 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});
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", "VELOCITY", "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<VELOCITY>(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, VELOCITY, 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
294 BOOST_AUTO_TEST_CASE(Upwind_gradient_3D_vector_velocity_test)
295 {
296 {
297 const size_t grid_dim = 3;
298 const double box_lower = -1.0;
299 const double box_upper = 1.0;
300 Box<grid_dim, double> box({box_lower, box_lower, box_lower}, {box_upper, box_upper, box_upper});
303 typedef grid_dist_id<grid_dim, double, props> grid_in_type;
304
305 double mu = 0.5 * (box_upper - abs(box_lower));
306 double sigma = 0.3 * (box_upper - box_lower);
307
308 size_t N = 32;
309// for(size_t N = 32; N <=256; N *=2)
310 {
311 const size_t sz[grid_dim] = {N, N, N};
312 grid_in_type g_dist(sz, box, ghost);
313 g_dist.setPropNames({"F_GAUSSIAN", "Velocity", "dF_GAUSSIAN", "dF_UPWIND", "ERROR"});
314
315 auto gdom = g_dist.getDomainGhostIterator();
316 while (gdom.isNext())
317 {
318 auto key = gdom.get();
319 Point<grid_dim, double> p = g_dist.getPos(key);
320 // Initialize grid and ghost with gaussian function
321 g_dist.getProp<F_GAUSSIAN>(key) = gaussian(p, mu, sigma);
322 ++gdom;
323 }
324
325 auto dom = g_dist.getDomainIterator();
326 while (dom.isNext())
327 {
328 auto key = dom.get();
329 Point<grid_dim, double> p = g_dist.getPos(key);
330
331 for (int d = 0; d < grid_dim; d++)
332 {
333 // Analytical gradient
334 g_dist.getProp<dF_GAUSSIAN>(key)[d] = hermite_polynomial(p.get(d), sigma, 1) * g_dist.getProp<F_GAUSSIAN>(key);
335 // Initialize the velocity vector field
336 g_dist.getProp<VELOCITY>(key)[d] = g_dist.getProp<dF_GAUSSIAN>(key)[d];
337 }
338 ++dom;
339 }
340
341 for (int convergence_order = 1; convergence_order <=5; convergence_order +=2 )
342 {
343 // Upwind gradient with specific convergence order, not becoming one-sided at the boundary
344 get_upwind_gradient<F_GAUSSIAN, VELOCITY, dF_UPWIND>(g_dist, convergence_order, false);
345 // Get the error between analytical and numerical solution
346 get_relative_error<dF_UPWIND, dF_GAUSSIAN, ERROR>(g_dist);
347
348 LNorms<double> lNorms;
349 lNorms.get_l_norms_grid<ERROR>(g_dist);
350
351 if(N==32)
352 {
353 switch(convergence_order)
354 {
355 case 1:
356 BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.38954184748194 + EPSILON, "Checking L2-norm upwind"
357 " gradient with vector "
358 "velocity order " +
359 std::to_string(convergence_order));
360 break;
361 case 3:
362 BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.08667855716144 + EPSILON, "Checking L2-norm upwind"
363 " gradient with "
364 "vector velocity order " +
365 std::to_string(convergence_order));
366 break;
367 case 5:
368 BOOST_CHECK_MESSAGE(lNorms.l2 <= 0.02285996528578 + EPSILON, "Checking L2-norm upwind"
369 " gradient with "
370 "vector velocity order " +
371 std::to_string(convergence_order));
372
373 break;
374 default:
375 std::cout << "Checking only implemented for convergence order 1, 3 and 5." << std::endl;
376 }
377 }
378 }
379
380 }
381 }
382 }
383BOOST_AUTO_TEST_SUITE_END()
Header file containing small mathematical help-functions that are needed e.g. for the Sussman redista...
int sgn(T val)
Gets the sign of a variable.
Header file containing functions for computing the error and the L_2 / L_infinity norm.
Approximating upwind gradients on a grid with the following options for the order of accuracy: 1,...
This class represent an N-dimensional box.
Definition Box.hpp:61
Class for computing the l2/l_infinity norm for distributed grids and vectors based on given errors.
Definition LNorms.hpp:105
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
This class implement the point shape in an N-dimensional space.
Definition Point.hpp:28
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition Point.hpp:172
This is a distributed grid.
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...