1/*
2 * eq_unit_test_3d.hpp
3 *
4 * Created on: Jan 4, 2016
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_3D_HPP_
9#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_3D_HPP_
10
11#define BOOST_TEST_DYN_LINK
12#include <boost/test/unit_test.hpp>
13
14#include "config.h"
15#include "Laplacian.hpp"
16#include "FiniteDifference/eq.hpp"
17#include "FiniteDifference/sum.hpp"
18#include "FiniteDifference/mul.hpp"
19#include "Grid/grid_dist_id.hpp"
20#include "Decomposition/CartDecomposition.hpp"
21#include "Vector/Vector.hpp"
22#include "Solvers/umfpack_solver.hpp"
23#include "data_type/aggregate.hpp"
24#include "Solvers/petsc_solver.hpp"
25#include "FiniteDifference/FDScheme.hpp"
26
27BOOST_AUTO_TEST_SUITE( eq_test_suite_3d )
28
29//! Specify the general caratteristic of system to solve
30struct lid_nn_3d_eigen
31{
32 //! dimensionaly of the equation ( 3D problem ...)
33 static const unsigned int dims = 3;
34 //! number of fields in the system
35 static const unsigned int nvar = 4;
36
37 //! boundary at X and Y
38 static const bool boundary[];
39
40 //! type of space float, double, ...
41 typedef float stype;
42
43 //! type of base grid
44 typedef grid_dist_id<3,float,aggregate<float[3],float>,CartDecomposition<3,float>> b_grid;
45
46 //! type of SparseMatrix for the linear solver
47 typedef SparseMatrix<double,int,EIGEN_BASE> SparseMatrix_type;
48
49 //! type of Vector for the linear solver
50 typedef Vector<double> Vector_type;
51
52 //! Define the underline grid is staggered
53 static const int grid_type = STAGGERED_GRID;
54};
55
56struct lid_nn_3d_petsc
57{
58 //! dimensionaly of the equation ( 3D problem ...)
59 static const unsigned int dims = 3;
60 //! number of fields in the system
61 static const unsigned int nvar = 4;
62
63 //! boundary at X and Y
64 static const bool boundary[];
65
66 //! type of space float, double, ...
67 typedef float stype;
68
69 //! type of base grid
70 typedef grid_dist_id<3,float,aggregate<float[3],float>,CartDecomposition<3,float>> b_grid;
71
72 //! type of SparseMatrix for the linear solver
73 typedef SparseMatrix<double,int,PETSC_BASE> SparseMatrix_type;
74
75 //! type of Vector for the linear solver
76 typedef Vector<double,PETSC_BASE> Vector_type;
77
78 //! Define the underline grid is staggered
79 static const int grid_type = STAGGERED_GRID;
80};
81
82//typedef lid_nn_3d_eigen lid_nn_3d;
83
84const bool lid_nn_3d_eigen::boundary[] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
85const bool lid_nn_3d_petsc::boundary[] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
86
87// Constant Field
88struct eta
89{
90 //! define that eta is a constant field
91 typedef void const_field;
92
93 //! therutn the value of the constant
94 static float val() {return 1.0;}
95};
96
97template<typename solver_type,typename lid_nn_3d> void lid_driven_cavity_3d()
98{
99 #include "Equations/stoke_flow_eq_3d.hpp"
100
101 Vcluster<> & v_cl = create_vcluster();
102
103 if (v_cl.getProcessingUnits() > 3)
104 return;
105
106 // Domain
107 Box<3,float> domain({0.0,0.0,0.0},{3.0,1.0,1.0});
108
109 // Ghost
110 Ghost<3,float> g(0.01);
111
112 long int sz[] = {36,12,12};
113 size_t szu[3];
114 szu[0] = (size_t)sz[0];
115 szu[1] = (size_t)sz[1];
116 szu[2] = (size_t)sz[2];
117
118 Padding<3> pd({1,1,1},{0,0,0});
119
120 // velocity in the grid is the property 0, pressure is the property 1
121 constexpr int velocity = 0;
122 constexpr int pressure = 1;
123
124 // Initialize openfpm
125 grid_dist_id<3,float,aggregate<float[3],float>,CartDecomposition<3,float>> g_dist(szu,domain,g);
126
127 // Ghost stencil
128 Ghost<3,long int> stencil_max(1);
129
130 // Distributed grid
131 FDScheme<lid_nn_3d> fd(pd,stencil_max,domain,g_dist);
132
133 // start and end of the bulk
134
135 fd.impose(ic_eq(),0.0, EQ_4, {0,0,0},{sz[0]-2,sz[1]-2,sz[2]-2},true);
136 fd.impose(Prs(), 0.0, EQ_4, {0,0,0},{0,0,0});
137 fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2,sz[2]-2});
138 fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2,sz[2]-2});
139 fd.impose(vz_eq(),0.0, EQ_3, {0,0,1},{sz[0]-2,sz[1]-2,sz[2]-2});
140
141 // v_x
142 // R L
143 fd.impose(v_x(),0.0, EQ_1, {0,0,0}, {0,sz[1]-2,sz[2]-2});
144 fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0,0},{sz[0]-1,sz[1]-2,sz[2]-2});
145
146 // T B
147 fd.impose(avg_y_vx_f(),0.0, EQ_1, {0,-1,0}, {sz[0]-1,-1,sz[2]-2});
148 fd.impose(avg_y_vx(),0.0, EQ_1, {0,sz[1]-1,0},{sz[0]-1,sz[1]-1,sz[2]-2});
149
150 // A F
151 fd.impose(avg_z_vx_f(),0.0, EQ_1, {0,-1,-1}, {sz[0]-1,sz[1]-1,-1});
152 fd.impose(avg_z_vx(),0.0, EQ_1, {0,-1,sz[2]-1},{sz[0]-1,sz[1]-1,sz[2]-1});
153
154 // v_y
155 // R L
156 fd.impose(avg_x_vy_f(),0.0, EQ_2, {-1,0,0}, {-1,sz[1]-1,sz[2]-2});
157 fd.impose(avg_x_vy(),1.0, EQ_2, {sz[0]-1,0,0},{sz[0]-1,sz[1]-1,sz[2]-2});
158
159 // T B
160 fd.impose(v_y(), 0.0, EQ_2, {0,0,0}, {sz[0]-2,0,sz[2]-2});
161 fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1,0},{sz[0]-2,sz[1]-1,sz[2]-2});
162
163 // F A
164 fd.impose(avg_z_vy(),0.0, EQ_2, {-1,0,sz[2]-1}, {sz[0]-1,sz[1]-1,sz[2]-1});
165 fd.impose(avg_z_vy_f(),0.0, EQ_2, {-1,0,-1}, {sz[0]-1,sz[1]-1,-1});
166
167 // v_z
168 // R L
169 fd.impose(avg_x_vz_f(),0.0, EQ_3, {-1,0,0}, {-1,sz[1]-2,sz[2]-1});
170 fd.impose(avg_x_vz(),1.0, EQ_3, {sz[0]-1,0,0},{sz[0]-1,sz[1]-2,sz[2]-1});
171
172 // T B
173 fd.impose(avg_y_vz(),0.0, EQ_3, {-1,sz[1]-1,0},{sz[0]-1,sz[1]-1,sz[2]-1});
174 fd.impose(avg_y_vz_f(),0.0, EQ_3, {-1,-1,0}, {sz[0]-1,-1,sz[2]-1});
175
176 // F A
177 fd.impose(v_z(),0.0, EQ_3, {0,0,0}, {sz[0]-2,sz[1]-2,0});
178 fd.impose(v_z(),0.0, EQ_3, {0,0,sz[2]-1},{sz[0]-2,sz[1]-2,sz[2]-1});
179
180 // Padding pressure
181
182 // L R
183 fd.impose(Prs(), 0.0, EQ_4, {-1,-1,-1},{-1,sz[1]-1,sz[2]-1});
184 fd.impose(Prs(), 0.0, EQ_4, {sz[0]-1,-1,-1},{sz[0]-1,sz[1]-1,sz[2]-1});
185
186 // T B
187 fd.impose(Prs(), 0.0, EQ_4, {0,sz[1]-1,-1}, {sz[0]-2,sz[1]-1,sz[2]-1});
188 fd.impose(Prs(), 0.0, EQ_4, {0,-1 ,-1}, {sz[0]-2,-1, sz[2]-1});
189
190 // F A
191 fd.impose(Prs(), 0.0, EQ_4, {0,0,sz[2]-1}, {sz[0]-2,sz[1]-2,sz[2]-1});
192 fd.impose(Prs(), 0.0, EQ_4, {0,0,-1}, {sz[0]-2,sz[1]-2,-1});
193
194 // Impose v_x v_y v_z padding
195 fd.impose(v_x(), 0.0, EQ_1, {-1,-1,-1},{-1,sz[1]-1,sz[2]-1});
196 fd.impose(v_y(), 0.0, EQ_2, {-1,-1,-1},{sz[0]-1,-1,sz[2]-1});
197 fd.impose(v_z(), 0.0, EQ_3, {-1,-1,-1},{sz[0]-1,sz[1]-1,-1});
198
199 solver_type solver;
200 auto x_ = solver.try_solve(fd.getA(),fd.getB());
201
202 // Bring the solution to grid
203 fd.template copy<velocity,pressure>(x_,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
204
205 std::string s = std::string(demangle(typeid(solver_type).name()));
206 s += "_";
207
208 g_dist.write(s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid");
209
210#if !(defined(SE_CLASS3) || defined(TEST_COVERAGE_MODE))
211
212 // Initialize openfpm
213 grid_dist_id<3,float,aggregate<float[3],float>,CartDecomposition<3,float>> g_dist2(g_dist.getDecomposition(),szu,g);
214 g_dist2.load("test/lid_driven_cavity_3d_reference.hdf5");
215
216 auto it2 = g_dist2.getDomainIterator();
217
218 bool test = true;
219 while (it2.isNext())
220 {
221 auto p = it2.get();
222
223 test &= fabs(g_dist2.template getProp<velocity>(p)[0] - g_dist.template getProp<velocity>(p)[0]) < 3.5e-5;
224 test &= fabs(g_dist2.template getProp<velocity>(p)[1] - g_dist.template getProp<velocity>(p)[1]) < 3.5e-5;
225 test &= fabs(g_dist2.template getProp<velocity>(p)[2] - g_dist.template getProp<velocity>(p)[2]) < 3.5e-5;
226
227 test &= fabs(g_dist2.template getProp<pressure>(p) - g_dist.template getProp<pressure>(p)) < 3.0e-4;
228
229 if (test == false)
230 {
231 std::cout << g_dist2.template getProp<velocity>(p)[0] << " " << g_dist.template getProp<velocity>(p)[0] << std::endl;
232 std::cout << g_dist2.template getProp<velocity>(p)[1] << " " << g_dist.template getProp<velocity>(p)[1] << std::endl;
233 std::cout << g_dist2.template getProp<velocity>(p)[2] << " " << g_dist.template getProp<velocity>(p)[2] << std::endl;
234
235 std::cout << g_dist2.template getProp<pressure>(p) << " " << g_dist.template getProp<pressure>(p) << std::endl;
236
237 break;
238 }
239
240 ++it2;
241 }
242
243 BOOST_REQUIRE_EQUAL(test,true);
244
245#endif
246
247}
248
249// Lid driven cavity, uncompressible fluid
250
251BOOST_AUTO_TEST_CASE(lid_driven_cavity)
252{
253#if defined(HAVE_EIGEN) && defined(HAVE_SUITESPARSE)
254 lid_driven_cavity_3d<umfpack_solver<double>,lid_nn_3d_eigen>();
255#endif
256#ifdef HAVE_PETSC
257 lid_driven_cavity_3d<petsc_solver<double>,lid_nn_3d_petsc>();
258#endif
259}
260
261BOOST_AUTO_TEST_SUITE_END()
262
263
264#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_3D_HPP_ */
265