OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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
14BOOST_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 }
206BOOST_AUTO_TEST_SUITE_END()
HelpFunctions.hpp
Header file containing small mathematical help-functions that are needed e.g. for the Sussman redista...
LNorms.hpp
Header file containing functions for computing the error and the L_2 / L_infinity norm.
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
Field
Definition eq.hpp:83
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
Error
Out-of-bound policy kill the program.
Definition vector_dist_ofb.hpp:70
aggregate
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition aggregate.hpp:215