OpenFPM_pdata  1.1.0
Project that contain the implementation of distributed structures
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Friends Pages
eq_unit_test.hpp
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 #include "Laplacian.hpp"
12 #include "FiniteDifference/eq.hpp"
13 #include "FiniteDifference/sum.hpp"
14 #include "FiniteDifference/mul.hpp"
15 #include "Grid/grid_dist_id.hpp"
16 #include "Decomposition/CartDecomposition.hpp"
17 #include "Vector/Vector.hpp"
18 #include "Solvers/umfpack_solver.hpp"
19 #include "data_type/aggregate.hpp"
20 #include "FiniteDifference/FDScheme.hpp"
21 
22 BOOST_AUTO_TEST_SUITE( eq_test_suite )
23 
24 
26 struct lid_nn
27 {
28  // dimensionaly of the equation (2D problem 3D problem ...)
29  static const unsigned int dims = 2;
30 
31  // number of fields in the system v_x, v_y, P so a total of 3
32  static const unsigned int nvar = 3;
33 
34  // boundary conditions PERIODIC OR NON_PERIODIC
35  static const bool boundary[];
36 
37  // type of space float, double, ...
38  typedef float stype;
39 
40  // type of base grid, it is the distributed grid that will store the result
41  // Note the first property is a 2D vector (velocity), the second is a scalar (Pressure)
43 
44  // type of SparseMatrix, for the linear system, this parameter is bounded by the solver
45  // that you are using, in case of umfpack it is the only possible choice
47 
48  // type of Vector for the linear system, this parameter is bounded by the solver
49  // that you are using, in case of umfpack it is the only possible choice
51 
52  // Define that the underline grid where we discretize the system of equation is staggered
53  static const int grid_type = STAGGERED_GRID;
54 };
55 
56 const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC};
57 
59 
61 
62 // Constant Field
63 struct eta
64 {
65  typedef void const_field;
66 
67  static float val() {return 1.0;}
68 };
69 
70 // Convenient constants
71 constexpr unsigned int v[] = {0,1};
72 constexpr unsigned int P = 2;
73 constexpr unsigned int ic = 2;
74 
75 // Create field that we have v_x, v_y, P
76 typedef Field<v[x],lid_nn> v_x;
77 typedef Field<v[y],lid_nn> v_y;
78 typedef Field<P,lid_nn> Prs;
79 
80 // Eq1 V_x
81 
83 typedef D<x,Prs,lid_nn> p_x;
84 typedef minus<p_x,lid_nn> m_p_x;
86 
87 // Eq2 V_y
88 
90 typedef D<y,Prs,lid_nn> p_y;
91 typedef minus<p_y,lid_nn> m_p_y;
93 
94 // Eq3 Incompressibility
95 
99 
100 
101 // Equation for boundary conditions
102 
103 /* Consider the staggered cell
104  *
105  \verbatim
106 
107  +--$--+
108  | |
109  # * #
110  | |
111  0--$--+
112 
113  # = velocity(y)
114  $ = velocity(x)
115  * = pressure
116 
117  \endverbatim
118  *
119  *
120  * If we want to impose v_y = 0 on 0 we have to interpolate between # of this cell
121  * and # of the previous cell on y, (Average) or Avg operator
122  *
123  */
124 
125 // Directional Avg
126 typedef Avg<x,v_y,lid_nn> avg_vy;
127 typedef Avg<y,v_x,lid_nn> avg_vx;
128 
131 
132 #define EQ_1 0
133 #define EQ_2 1
134 #define EQ_3 2
135 
137 
138 template<typename solver_type,typename lid_nn> void lid_driven_cavity_2d()
139 {
140  Vcluster & v_cl = create_vcluster();
141 
142  if (v_cl.getProcessingUnits() > 3)
143  return;
144 
146 
147  // velocity in the grid is the property 0, pressure is the property 1
148  constexpr int velocity = 0;
149  constexpr int pressure = 1;
150 
151  // Domain, a rectangle
152  Box<2,float> domain({0.0,0.0},{3.0,1.0});
153 
154  // Ghost (Not important in this case but required)
155  Ghost<2,float> g(0.01);
156 
157  // Grid points on x=256 and y=64
158  long int sz[] = {256,64};
159  size_t szu[2];
160  szu[0] = (size_t)sz[0];
161  szu[1] = (size_t)sz[1];
162 
163  // We need one more point on the left and down part of the domain
164  // This is given by the boundary conditions that we impose, the
165  // reason is mathematical in order to have a well defined system
166  // and cannot be discussed here
167  Padding<2> pd({1,1},{0,0});
168 
169  // Distributed grid that store the solution
171 
172  // It is the maximum extension of the stencil
173  Ghost<2,long int> stencil_max(1);
174 
175  // Finite difference scheme
176  FDScheme<lid_nn> fd(pd, stencil_max, domain,g_dist);
177 
178  // Here we impose the equation, we start from the incompressibility Eq imposed in the bulk with the
179  // exception of the first point {0,0} and than we set P = 0 in {0,0}, why we are doing this is again
180  // mathematical to have a well defined system, an intuitive explanation is that P and P + c are both
181  // solution for the incompressibility equation, this produce an ill-posed problem to make it well posed
182  // we set one point in this case {0,0} the pressure to a fixed constant for convenience P = 0
183  fd.impose(ic_eq(),0.0, EQ_3, {0,0},{sz[0]-2,sz[1]-2},true);
184  fd.impose(Prs(), 0.0, EQ_3, {0,0},{0,0});
185 
186  // Here we impose the Eq1 and Eq2
187  fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2});
188  fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2});
189 
190  // v_x and v_y
191  // Imposing B1
192  fd.impose(v_x(),0.0, EQ_1, {0,0},{0,sz[1]-2});
193  fd.impose(avg_vy_f(),0.0, EQ_2 , {-1,0},{-1,sz[1]-1});
194  // Imposing B2
195  fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
196  fd.impose(avg_vy(),1.0, EQ_2, {sz[0]-1,0},{sz[0]-1,sz[1]-1});
197 
198  // Imposing B3
199  fd.impose(avg_vx_f(),0.0, EQ_1, {0,-1},{sz[0]-1,-1});
200  fd.impose(v_y(), 0.0, EQ_2, {0,0},{sz[0]-2,0});
201  // Imposing B4
202  fd.impose(avg_vx(),0.0, EQ_1, {0,sz[1]-1},{sz[0]-1,sz[1]-1});
203  fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1},{sz[0]-2,sz[1]-1});
204 
205  // When we pad the grid, there are points of the grid that are not
206  // touched by the previous condition. Mathematically this lead
207  // to have too many variables for the conditions that we are imposing.
208  // Here we are imposing variables that we do not touch to zero
209  //
210 
211  // Padding pressure
212  fd.impose(Prs(), 0.0, EQ_3, {-1,-1},{sz[0]-1,-1});
213  fd.impose(Prs(), 0.0, EQ_3, {-1,sz[1]-1},{sz[0]-1,sz[1]-1});
214  fd.impose(Prs(), 0.0, EQ_3, {-1,0},{-1,sz[1]-2});
215  fd.impose(Prs(), 0.0, EQ_3, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
216 
217  // Impose v_x Padding Impose v_y padding
218  fd.impose(v_x(), 0.0, EQ_1, {-1,-1},{-1,sz[1]-1});
219  fd.impose(v_y(), 0.0, EQ_2, {-1,-1},{sz[0]-1,-1});
220 
221  solver_type solver;
222  auto x = solver.solve(fd.getA(),fd.getB());
223 
225 
227 
228  fd.template copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1},g_dist);
229 
230  std::string s = std::string(demangle(typeid(solver_type).name()));
231  s += "_";
232 
234 
235  g_dist.write(s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid");
236 
237 #ifdef HAVE_OSX
238 
239  std::string file1 = std::string("test/") + s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_osx.vtk";
240  std::string file2 = s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
241 
242 #else
243 
244  #if __GNUC__ == 6
245 
246  std::string file1 = std::string("test/") + s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC6.vtk";
247  std::string file2 = s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
248 
249  #elif __GNUC__ == 5
250 
251  std::string file1 = std::string("test/") + s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC5.vtk";
252  std::string file2 = s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
253 
254  #else
255 
256  std::string file1 = std::string("test/") + s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC4.vtk";
257  std::string file2 = s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
258 
259  #endif
260 
261 
262 #endif
263 
264  std::cout << "File1: " << file1 << std::endl;
265  std::cout << "File2: " << file2 << std::endl;
266 
267 #ifndef SE_CLASS3
268 
269  // Check that match
270  bool test = compare(file1,file2);
271  BOOST_REQUIRE_EQUAL(test,true);
272 
273 #endif
274 
275 }
276 
277 // Lid driven cavity, incompressible fluid
278 
279 BOOST_AUTO_TEST_CASE(lid_driven_cavity)
280 {
281  lid_driven_cavity_2d<umfpack_solver<double>,lid_nn>();
282 }
283 
284 BOOST_AUTO_TEST_SUITE_END()
285 
286 #endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_HPP_ */
Derivative second order on h (spacing)
Definition: Derivative.hpp:28
[Definition of the system]
It ancapsulate the minus operation.
Definition: sum.hpp:141
size_t getProcessUnitID()
Get the process unit id.
Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support...
Definition: Vector.hpp:39
Definition: eq.hpp:81
Class that contain Padding information on each direction positive and Negative direction.
Definition: Ghost.hpp:117
Definition: Ghost.hpp:39
Implementation of VCluster class.
Definition: VCluster.hpp:36
It model an expression expr1 * expr2.
Definition: mul.hpp:119
This class decompose a space into sub-sub-domains and distribute them across processors.
Sparse Matrix implementation.
This is a distributed grid.
Average.
Definition: Average.hpp:31
[Definition of the system]
This class represent an N-dimensional box.
Definition: Box.hpp:56
It model an expression expr1 + ... exprn.
Definition: sum.hpp:92
Finite Differences.
Definition: FDScheme.hpp:124
Test structure used for several test.
Definition: Point_test.hpp:105
size_t getProcessingUnits()
Get the total number of processors.