1/*
2 * eq_unit_test.hpp
3 *
4 * Created on: Oct 13, 2015
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_HPP_
9#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_HPP_
10
11#define BOOST_TEST_DYN_LINK
12#include <boost/test/unit_test.hpp>
13
14#include "Laplacian.hpp"
15#include "FiniteDifference/eq.hpp"
16#include "FiniteDifference/sum.hpp"
17#include "FiniteDifference/mul.hpp"
18#include "Grid/grid_dist_id.hpp"
19#include "Decomposition/CartDecomposition.hpp"
20#include "Vector/Vector.hpp"
21#include "Solvers/umfpack_solver.hpp"
22#include "data_type/aggregate.hpp"
23#include "FiniteDifference/FDScheme.hpp"
24
25constexpr unsigned int x = 0;
26constexpr unsigned int y = 1;
27constexpr unsigned int z = 2;
28constexpr unsigned int V = 0;
29
30BOOST_AUTO_TEST_SUITE( eq_test_suite )
31
32//! [Definition of the system]
33
34struct lid_nn
35{
36 // dimensionaly of the equation (2D problem 3D problem ...)
37 static const unsigned int dims = 2;
38
39 // number of fields in the system v_x, v_y, P so a total of 3
40 static const unsigned int nvar = 3;
41
42 // boundary conditions PERIODIC OR NON_PERIODIC
43 static const bool boundary[];
44
45 // type of space float, double, ...
46 typedef float stype;
47
48 // type of base grid, it is the distributed grid that will store the result
49 // Note the first property is a 2D vector (velocity), the second is a scalar (Pressure)
50 typedef grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> b_grid;
51
52 // type of SparseMatrix, for the linear system, this parameter is bounded by the solver
53 // that you are using, in case of umfpack it is the only possible choice
54 typedef SparseMatrix<double,int> SparseMatrix_type;
55
56 // type of Vector for the linear system, this parameter is bounded by the solver
57 // that you are using, in case of umfpack it is the only possible choice
58 typedef Vector<double> Vector_type;
59
60 // Define that the underline grid where we discretize the system of equation is staggered
61 static const int grid_type = STAGGERED_GRID;
62};
63
64const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC};
65
66//! [Definition of the system]
67
68//! [Definition of the equation of the system in the bulk and at the boundary]
69
70// Constant Field
71struct eta
72{
73 typedef void const_field;
74
75 static float val() {return 1.0;}
76};
77
78// Convenient constants
79constexpr unsigned int v[] = {0,1};
80constexpr unsigned int P = 2;
81constexpr unsigned int ic = 2;
82
83// Create field that we have v_x, v_y, P
84typedef Field<v[x],lid_nn> v_x;
85typedef Field<v[y],lid_nn> v_y;
86typedef Field<P,lid_nn> Prs;
87
88// Eq1 V_x
89
90typedef mul<eta,Lap<v_x,lid_nn>,lid_nn> eta_lap_vx;
91typedef D<x,Prs,lid_nn> p_x;
92typedef minus<p_x,lid_nn> m_p_x;
93typedef sum<eta_lap_vx,m_p_x,lid_nn> vx_eq;
94
95// Eq2 V_y
96
97typedef mul<eta,Lap<v_y,lid_nn>,lid_nn> eta_lap_vy;
98typedef D<y,Prs,lid_nn> p_y;
99typedef minus<p_y,lid_nn> m_p_y;
100typedef sum<eta_lap_vy,m_p_y,lid_nn> vy_eq;
101
102// Eq3 Incompressibility
103
104typedef D<x,v_x,lid_nn,FORWARD> dx_vx;
105typedef D<y,v_y,lid_nn,FORWARD> dy_vy;
106typedef sum<dx_vx,dy_vy,lid_nn> ic_eq;
107
108
109// Equation for boundary conditions
110
111/* Consider the staggered cell
112 *
113 \verbatim
114
115 +--$--+
116 | |
117 # * #
118 | |
119 0--$--+
120
121 # = velocity(y)
122 $ = velocity(x)
123 * = pressure
124
125 \endverbatim
126 *
127 *
128 * If we want to impose v_y = 0 on 0 we have to interpolate between # of this cell
129 * and # of the previous cell on y, (Average) or Avg operator
130 *
131 */
132
133// Directional Avg
134typedef Avg<x,v_y,lid_nn> avg_vy;
135typedef Avg<y,v_x,lid_nn> avg_vx;
136
137typedef Avg<x,v_y,lid_nn,FORWARD> avg_vy_f;
138typedef Avg<y,v_x,lid_nn,FORWARD> avg_vx_f;
139
140#define EQ_1 0
141#define EQ_2 1
142#define EQ_3 2
143
144//! [Definition of the equation of the system in the bulk and at the boundary]
145
146template<typename solver_type,typename lid_nn> void lid_driven_cavity_2d()
147{
148 Vcluster<> & v_cl = create_vcluster();
149
150 if (v_cl.getProcessingUnits() > 3)
151 return;
152
153 //! [lid-driven cavity 2D]
154
155 // velocity in the grid is the property 0, pressure is the property 1
156 constexpr int velocity = 0;
157 constexpr int pressure = 1;
158
159 // Domain, a rectangle
160 Box<2,float> domain({0.0,0.0},{3.0,1.0});
161
162 // Ghost (Not important in this case but required)
163 Ghost<2,float> g(0.01);
164
165 // Grid points on x=256 and y=64
166 long int sz[] = {256,64};
167 size_t szu[2];
168 szu[0] = (size_t)sz[0];
169 szu[1] = (size_t)sz[1];
170
171 // We need one more point on the left and down part of the domain
172 // This is given by the boundary conditions that we impose, the
173 // reason is mathematical in order to have a well defined system
174 // and cannot be discussed here
175 Padding<2> pd({1,1},{0,0});
176
177 // Distributed grid that store the solution
178 grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> g_dist(szu,domain,g);
179
180 // It is the maximum extension of the stencil
181 Ghost<2,long int> stencil_max(1);
182
183 // Finite difference scheme
184 FDScheme<lid_nn> fd(pd, stencil_max, domain,g_dist);
185
186 // Here we impose the equation, we start from the incompressibility Eq imposed in the bulk with the
187 // exception of the first point {0,0} and than we set P = 0 in {0,0}, why we are doing this is again
188 // mathematical to have a well defined system, an intuitive explanation is that P and P + c are both
189 // solution for the incompressibility equation, this produce an ill-posed problem to make it well posed
190 // we set one point in this case {0,0} the pressure to a fixed constant for convenience P = 0
191 fd.impose(ic_eq(),0.0, EQ_3, {0,0},{sz[0]-2,sz[1]-2},true);
192 fd.impose(Prs(), 0.0, EQ_3, {0,0},{0,0});
193
194 // Here we impose the Eq1 and Eq2
195 fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2});
196 fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2});
197
198 // v_x and v_y
199 // Imposing B1
200 fd.impose(v_x(),0.0, EQ_1, {0,0},{0,sz[1]-2});
201 fd.impose(avg_vy_f(),0.0, EQ_2 , {-1,0},{-1,sz[1]-1});
202 // Imposing B2
203 fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
204 fd.impose(avg_vy(),1.0, EQ_2, {sz[0]-1,0},{sz[0]-1,sz[1]-1});
205
206 // Imposing B3
207 fd.impose(avg_vx_f(),0.0, EQ_1, {0,-1},{sz[0]-1,-1});
208 fd.impose(v_y(), 0.0, EQ_2, {0,0},{sz[0]-2,0});
209 // Imposing B4
210 fd.impose(avg_vx(),0.0, EQ_1, {0,sz[1]-1},{sz[0]-1,sz[1]-1});
211 fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1},{sz[0]-2,sz[1]-1});
212
213 // When we pad the grid, there are points of the grid that are not
214 // touched by the previous condition. Mathematically this lead
215 // to have too many variables for the conditions that we are imposing.
216 // Here we are imposing variables that we do not touch to zero
217 //
218
219 // Padding pressure
220 fd.impose(Prs(), 0.0, EQ_3, {-1,-1},{sz[0]-1,-1});
221 fd.impose(Prs(), 0.0, EQ_3, {-1,sz[1]-1},{sz[0]-1,sz[1]-1});
222 fd.impose(Prs(), 0.0, EQ_3, {-1,0},{-1,sz[1]-2});
223 fd.impose(Prs(), 0.0, EQ_3, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
224
225 // Impose v_x Padding Impose v_y padding
226 fd.impose(v_x(), 0.0, EQ_1, {-1,-1},{-1,sz[1]-1});
227 fd.impose(v_y(), 0.0, EQ_2, {-1,-1},{sz[0]-1,-1});
228
229 solver_type solver;
230 auto x = solver.solve(fd.getA(),fd.getB());
231
232 //! [lid-driven cavity 2D]
233
234 //! [Copy the solution to grid]
235
236 fd.template copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1},g_dist);
237
238 std::string s = std::string(demangle(typeid(solver_type).name()));
239 s += "_";
240
241 //! [Copy the solution to grid]
242
243 g_dist.write(s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid");
244
245#if !(defined(SE_CLASS3) || defined(TEST_COVERAGE_MODE))
246
247 // Initialize openfpm
248 grid_dist_id<2,float,aggregate<float[2],float>> g_dist2(g_dist.getDecomposition(),szu,g);
249 g_dist2.load("test/lid_driven_cavity_reference.hdf5");
250
251 auto it2 = g_dist2.getDomainIterator();
252
253 bool test = true;
254 while (it2.isNext())
255 {
256 auto p = it2.get();
257
258 test &= fabs(g_dist2.template getProp<velocity>(p)[0] - g_dist.template getProp<velocity>(p)[0]) < 3.5e-5;
259 test &= fabs(g_dist2.template getProp<velocity>(p)[1] - g_dist.template getProp<velocity>(p)[1]) < 3.5e-5;
260
261 test &= fabs(g_dist2.template getProp<pressure>(p) - g_dist.template getProp<pressure>(p)) < 3.0e-4;
262
263 if (test == false)
264 {
265 std::cout << g_dist2.template getProp<velocity>(p)[0] << " " << g_dist.template getProp<velocity>(p)[0] << std::endl;
266 std::cout << g_dist2.template getProp<velocity>(p)[1] << " " << g_dist.template getProp<velocity>(p)[1] << std::endl;
267
268 std::cout << g_dist2.template getProp<pressure>(p) << " " << g_dist.template getProp<pressure>(p) << std::endl;
269
270 break;
271 }
272
273 ++it2;
274 }
275
276 BOOST_REQUIRE_EQUAL(test,true);
277
278#endif
279}
280
281// Lid driven cavity, incompressible fluid
282
283BOOST_AUTO_TEST_CASE(lid_driven_cavity)
284{
285#if defined(HAVE_EIGEN) && defined(HAVE_SUITESPARSE)
286 lid_driven_cavity_2d<umfpack_solver<double>,lid_nn>();
287#endif
288}
289
290BOOST_AUTO_TEST_SUITE_END()
291
292#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_HPP_ */
293