OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
convergence_test.cpp
1//
2// Created by jstark on 27.05.21.
3//
4#define BOOST_TEST_DYN_LINK
5//#define BOOST_TEST_MAIN // in only one cpp file
6#include <boost/test/unit_test.hpp>
7
8// Include redistancing files
12// Include header files for testing
13#include "Draw/DrawSphere.hpp"
14#include "l_norms/LNorms.hpp"
16
17
18
19BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
20 typedef double phi_type;
21 typedef double space_type;
22 const phi_type EPSILON = std::numeric_limits<phi_type>::epsilon();
23
24 space_type get_min_required_time_step(size_t Nmax, space_type b_low, space_type b_up, int dims)
25 {
26 space_type dx = (b_up - b_low) / (space_type) Nmax;
27 space_type sum = 0;
28 for (int d = 0; d < dims; d++)
29 {
30 sum += 1 / dx;
31 }
32 return 0.1 / sum;
33 }
34
35 BOOST_AUTO_TEST_CASE(RedistancingSussman_convergence_test_1D)
36 {
37 // CFL dt in 1D for N=64 and c=1 is 0.000503905
38
39 const size_t grid_dim = 1;
40 // some indices
41 const size_t x = 0;
42
43 const size_t Phi_0_grid = 0;
44 const size_t SDF_sussman_grid = 1;
45 const size_t SDF_exact_grid = 2;
46 const size_t Error_grid = 3;
47
48 const space_type box_lower = -1.0;
49 const space_type box_upper = 1.0;
50 Box<grid_dim, space_type> box({box_lower}, {box_upper});
54
55// space_type dt = 0.000503905;
56 size_t Nmin = 32;
57 size_t Nmax = 4096;
58 space_type dt = get_min_required_time_step(Nmax, box_lower, box_upper, grid_dim);
59 for (size_t N = Nmin; N <= Nmax; N*=2)
60 {
61 const size_t sz[grid_dim] = {N};
62 grid_in_type g_dist(sz, box, ghost);
63 g_dist.setPropNames({"Phi_0", "SDF_sussman", "SDF_exact", "Relative error"});
64
65 // Now we initialize the grid with the indicator function and the analytical solution
66 auto dom = g_dist.getDomainIterator();
67 while(dom.isNext())
68 {
69 auto key = dom.get();
70 Point<grid_dim, space_type> coords = g_dist.getPos(key);
71 g_dist.template get<Phi_0_grid>(key) = sgn(coords.get(x));
72 g_dist.template get<SDF_exact_grid>(key) = coords.get(x);
73 ++dom;
74 }
75
76 {
77 Redist_options<phi_type> redist_options;
78 redist_options.min_iter = 1e3;
79 redist_options.max_iter = 1e16;
80
81 // set both convergence criteria to false s.t. termination only when max_iterations reached
82 redist_options.convTolChange.check = true; // define here which of the convergence criteria above should
83 // be used. If both are true, termination only occurs when both are fulfilled or when iter > max_iter
84 redist_options.convTolChange.value = EPSILON;
85
86 redist_options.convTolResidual.check = false; // (default: false)
87 redist_options.interval_check_convergence = 1000; // interval of #iterations at which
88 // convergence is checked (default: 100)
89 redist_options.width_NB_in_grid_points = 10; // width of narrow band in number of grid points.
90 // Must
91 // be at least 4, in order to have at least 2 grid points on each side of the interface. (default: 4)
92 redist_options.print_current_iterChangeResidual = true; // if true, prints out every current iteration + corresponding change from the previous iteration + residual from SDF (default: false)
93 redist_options.print_steadyState_iter = true; // if true, prints out the final iteration number when steady state was reached + final change + residual (default: true)
94
96 redist_options); // Instantiation of Sussman-redistancing class
97 std::cout << "CFL dt = " << redist_obj.get_time_step() << std::endl;
98 redist_obj.set_user_time_step(dt);
99 std::cout << "dt set to = " << dt << std::endl;
100 // Run the redistancing. in the <> brackets provide property-index where 1.) your initial Phi is stored and 2.) where the resulting SDF should be written to.
101 redist_obj.run_redistancing<Phi_0_grid, SDF_sussman_grid>();
102
103 // Compute the absolute error between analytical and numerical solution at each grid point
104 get_absolute_error<SDF_sussman_grid, SDF_exact_grid, Error_grid>(g_dist);
105// g_dist.write("grid_postRedistancing", FORMAT_BINARY);
107 // Get narrow band: Place particles on interface (narrow band width e.g. 4 grid points on each side of the
108 // interface)
109 size_t bc[grid_dim] = {NON_PERIODIC};
110 typedef aggregate<phi_type> props_nb;
112 Ghost<grid_dim, space_type> ghost_vd(0);
113 vd_type vd_narrow_band(0, box, bc, ghost_vd);
114 vd_narrow_band.setPropNames({"error"});
115 const size_t Error_vd = 0;
116 // Compute the L_2- and L_infinity-norm and save to file
117 NarrowBand<grid_in_type, phi_type> narrowBand_points(g_dist, redist_options.width_NB_in_grid_points ); // Instantiation of
118 // NarrowBand class
119 narrowBand_points.get_narrow_band_copy_specific_property<SDF_sussman_grid, Error_grid, Error_vd>(g_dist,
120 vd_narrow_band);
121// vd_narrow_band.write("vd_error_N" + std::to_string(N), FORMAT_BINARY);
122// vd_narrow_band.save("test_data/output/vd_nb8p_error" + std::to_string(N) + ".bin");
123 // Compute the L_2- and L_infinity-norm and save to file
124 LNorms<phi_type> lNorms_vd;
125 lNorms_vd.get_l_norms_vector<Error_vd>(vd_narrow_band);
126 lNorms_vd.write_to_file(N, 6, "l_norms_", "./");
127 }
128
129 }
130 }
131
132 BOOST_AUTO_TEST_CASE(RedistancingSussman_convergence_test_3D)
133 {
134 const size_t grid_dim = 3;
135 // some indices
136 const size_t x = 0;
137 const size_t y = 1;
138 const size_t z = 2;
139
140 const size_t Phi_0_grid = 0;
141 const size_t SDF_sussman_grid = 1;
142 const size_t SDF_exact_grid = 2;
143 const size_t Error_grid = 3;
144
145 for (size_t N=32; N <=128; N*=2)
146 {
147 const size_t sz[grid_dim] = {N, N, N};
148 const space_type radius = 1.0;
149 const space_type box_lower = 0.0;
150 const space_type box_upper = 4.0 * radius;
151 Box<grid_dim, space_type> box({box_lower, box_lower, box_lower}, {box_upper, box_upper, box_upper});
155 grid_in_type g_dist(sz, box, ghost);
156 g_dist.setPropNames({"Phi_0", "SDF_sussman", "SDF_exact", "Relative error"});
157
158 const space_type center[grid_dim] = {0.5*(box_upper+box_lower), 0.5*(box_upper+box_lower), 0.5*(box_upper+box_lower)};
159 init_grid_with_sphere<Phi_0_grid>(g_dist, radius, center[x], center[y], center[z]); // Initialize sphere onto grid
160
161 int order = 1;
162 {
163 Redist_options<phi_type> redist_options;
164 redist_options.min_iter = 1e3;
165 redist_options.max_iter = 1e3;
166
167 // set both convergence criteria to false s.t. termination only when max_iterations reached
168 redist_options.convTolChange.check = false; // define here which of the convergence criteria above should be used. If both are true, termination only occurs when both are fulfilled or when iter > max_iter
169 redist_options.convTolResidual.check = false; // (default: false)
170
171 redist_options.interval_check_convergence = 100; // interval of #iterations at which
172 // convergence is checked (default: 100)
173 redist_options.width_NB_in_grid_points = 4; // width of narrow band in number of
174 // grid points. Must be at least 4, in order to have at least 2 grid points on each side of the interface. (default: 4)
175 redist_options.print_current_iterChangeResidual = true; // if true, prints out every current iteration + corresponding change from the previous iteration + residual from SDF (default: false)
176 redist_options.print_steadyState_iter = true; // if true, prints out the final iteration number when steady state was reached + final change + residual (default: true)
177
178 RedistancingSussman<grid_in_type, phi_type> redist_obj(g_dist, redist_options); // Instantiation of
179 // Sussman-redistancing class
180 //redist_obj.set_user_time_step(dt);
181 //std::cout << "dt set to = " << dt << std::endl;
182 // Run the redistancing. in the <> brackets provide property-index where 1.) your initial Phi is stored and 2.) where the resulting SDF should be written to.
183 redist_obj.run_redistancing<Phi_0_grid, SDF_sussman_grid>();
184
185 // Compute exact signed distance function at each grid point
186 init_analytic_sdf_sphere<SDF_exact_grid>(g_dist, radius, center[x], center[y], center[z]);
187
188 // Compute the absolute error between analytical and numerical solution at each grid point
189 get_absolute_error<SDF_sussman_grid, SDF_exact_grid, Error_grid>(g_dist);
190
191
193 // Get narrow band: Place particles on interface (narrow band width e.g. 4 grid points on each side of the
194 // interface)
195 size_t bc[grid_dim] = {PERIODIC, PERIODIC, PERIODIC};
196 typedef aggregate<phi_type> props_nb;
198 Ghost<grid_dim, space_type> ghost_vd(0);
199 vd_type vd_narrow_band(0, box, bc, ghost_vd);
200 vd_narrow_band.setPropNames({"error"});
201 const size_t Error_vd = 0;
202 // Compute the L_2- and L_infinity-norm and save to file
203 NarrowBand<grid_in_type, phi_type> narrowBand_points(g_dist, redist_options.width_NB_in_grid_points); // Instantiation of NarrowBand class
204 narrowBand_points.get_narrow_band_copy_specific_property<SDF_sussman_grid, Error_grid, Error_vd>(g_dist,
205 vd_narrow_band);
206 vd_narrow_band.write("vd_nb4p_error_N" + std::to_string(N) + "_order" +
207 std::to_string(order), FORMAT_BINARY);
208// vd_narrow_band.save("test_data/output/vd_nb8p_error" + std::to_string(N) + ".bin");
209 // Compute the L_2- and L_infinity-norm and save to file
210 LNorms<phi_type> lNorms_vd;
211 lNorms_vd.get_l_norms_vector<Error_vd>(vd_narrow_band);
212 lNorms_vd.write_to_file(N, 6, "l_norms_vd_absError_4p_order" + std::to_string(order), "./");
213
214// 32,0.044692,0.066994
215// 64,0.010542,0.018396
216// 128,0.003301,0.009491
217
218 switch(N)
219 {
220 case 32:
221 BOOST_CHECK(lNorms_vd.l2 < 0.044693);
222 BOOST_CHECK(lNorms_vd.linf < 0.066995);
223 break;
224 case 64:
225 BOOST_CHECK(lNorms_vd.l2 < 0.010543);
226 BOOST_CHECK(lNorms_vd.linf < 0.018397);
227 break;
228 case 128:
229 BOOST_CHECK(lNorms_vd.l2 < 0.003302);
230 BOOST_CHECK(lNorms_vd.linf < 0.009492);
231 break;
232 }
233
234 }
235 }
236 }
237
238BOOST_AUTO_TEST_SUITE_END()
239
240
Header file containing functions that compute the analytical solution of the signed distance function...
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.
Class for getting the narrow band around the interface.
Header file containing functions for creating files and folders.
Class for reinitializing a level-set function into a signed distance function using Sussman redistanc...
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_vector(vectortype &vd)
Computes the L_2 and L_infinity norm on the basis of the precomputed error on a particle vector.
Definition LNorms.hpp:155
void write_to_file(const T col_0, const int precision, const std::string &filename, const std::string &path_output)
Writes the N (number of grid points on a square grid) and L-norms as strings to a csv-file.
Definition LNorms.hpp:226
Class for getting the narrow band around the interface.
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
Class for reinitializing a level-set function into a signed distance function using Sussman redistanc...
This is a distributed grid.
Distributed vector.
Structure to bundle options for redistancing.
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
It model an expression expr1 + ... exprn.
Definition sum.hpp:93