eq_unit_test.hpp 8.76 KB
Newer Older
incardon's avatar
incardon committed
1 2 3 4 5 6 7 8 9 10 11 12 13
/*
 * eq_unit_test.hpp
 *
 *  Created on: Oct 13, 2015
 *      Author: i-bird
 */

#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_HPP_

#include "Laplacian.hpp"
#include "FiniteDifference/eq.hpp"
#include "FiniteDifference/sum.hpp"
incardon's avatar
incardon committed
14
#include "FiniteDifference/mul.hpp"
incardon's avatar
incardon committed
15 16
#include "Grid/grid_dist_id.hpp"
#include "Decomposition/CartDecomposition.hpp"
17 18 19
#include "Vector/Vector.hpp"
#include "Solvers/umfpack_solver.hpp"
#include "data_type/aggregate.hpp"
20
#include "FiniteDifference/FDScheme.hpp"
incardon's avatar
incardon committed
21

22 23
BOOST_AUTO_TEST_SUITE( eq_test_suite )

24
//! [Definition of the system]
incardon's avatar
incardon committed
25 26 27 28 29

struct lid_nn
{
	// dimensionaly of the equation (2D problem 3D problem ...)
	static const unsigned int dims = 2;
30 31

	// number of fields in the system v_x, v_y, P so a total of 3
incardon's avatar
incardon committed
32 33
	static const unsigned int nvar = 3;

34
	// boundary conditions PERIODIC OR NON_PERIODIC
incardon's avatar
incardon committed
35 36 37 38
	static const bool boundary[];

	// type of space float, double, ...
	typedef float stype;
incardon's avatar
incardon committed
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)
42 43
	typedef grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> b_grid;

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
46 47
	typedef SparseMatrix<double,int> SparseMatrix_type;

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
50 51
	typedef Vector<double> Vector_type;

52
	// Define that the underline grid where we discretize the system of equation is staggered
53
	static const int grid_type = STAGGERED_GRID;
incardon's avatar
incardon committed
54 55 56 57
};

const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC};

58 59 60 61
//! [Definition of the system]

//! [Definition of the equation of the system in the bulk and at the boundary]

incardon's avatar
incardon committed
62 63 64
// Constant Field
struct eta
{
65 66 67
	typedef void const_field;

	static float val()	{return 1.0;}
incardon's avatar
incardon committed
68 69
};

70
// Convenient constants
incardon's avatar
incardon committed
71 72
constexpr unsigned int v[] = {0,1};
constexpr unsigned int P = 2;
73
constexpr unsigned int ic = 2;
incardon's avatar
incardon committed
74

75
// Create field that we have v_x, v_y, P
incardon's avatar
incardon committed
76 77 78 79 80 81 82 83
typedef Field<v[x],lid_nn> v_x;
typedef Field<v[y],lid_nn> v_y;
typedef Field<P,lid_nn> Prs;

// Eq1 V_x

typedef mul<eta,Lap<v_x,lid_nn>,lid_nn> eta_lap_vx;
typedef D<x,Prs,lid_nn> p_x;
84 85
typedef minus<p_x,lid_nn> m_p_x;
typedef sum<eta_lap_vx,m_p_x,lid_nn> vx_eq;
incardon's avatar
incardon committed
86 87 88 89

// Eq2 V_y

typedef mul<eta,Lap<v_y,lid_nn>,lid_nn> eta_lap_vy;
90 91 92
typedef D<y,Prs,lid_nn> p_y;
typedef minus<p_y,lid_nn> m_p_y;
typedef sum<eta_lap_vy,m_p_y,lid_nn> vy_eq;
incardon's avatar
incardon committed
93 94 95

// Eq3 Incompressibility

96 97
typedef D<x,v_x,lid_nn,FORWARD> dx_vx;
typedef D<y,v_y,lid_nn,FORWARD> dy_vy;
98 99
typedef sum<dx_vx,dy_vy,lid_nn> ic_eq;

incardon's avatar
incardon committed
100

101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
// Equation for boundary conditions

/* Consider the staggered cell
 *
 	 	 \verbatim

		+--$--+
		|     |
		#  *  #
		|     |
		0--$--+

	  # = velocity(y)
	  $ = velocity(x)
	  * = pressure

		\endverbatim
 *
 *
 * If we want to impose v_y = 0 on 0 we have to interpolate between # of this cell
 * and # of the previous cell on y, (Average) or Avg operator
 *
 */

125 126 127 128
// Directional Avg
typedef Avg<x,v_y,lid_nn> avg_vy;
typedef Avg<y,v_x,lid_nn> avg_vx;

129 130 131 132 133 134 135
typedef Avg<x,v_y,lid_nn,FORWARD> avg_vy_f;
typedef Avg<y,v_x,lid_nn,FORWARD> avg_vx_f;

#define EQ_1 0
#define EQ_2 1
#define EQ_3 2

136 137
//! [Definition of the equation of the system in the bulk and at the boundary]

Pietro Incardona's avatar
Pietro Incardona committed
138
template<typename solver_type,typename lid_nn> void lid_driven_cavity_2d()
incardon's avatar
incardon committed
139
{
140
	Vcluster & v_cl = create_vcluster();
141

142 143 144
	if (v_cl.getProcessingUnits() > 3)
		return;

145 146
	//! [lid-driven cavity 2D]

147 148 149 150
	// velocity in the grid is the property 0, pressure is the property 1
	constexpr int velocity = 0;
	constexpr int pressure = 1;

151
	// Domain, a rectangle
152
	Box<2,float> domain({0.0,0.0},{3.0,1.0});
incardon's avatar
incardon committed
153

154
	// Ghost (Not important in this case but required)
incardon's avatar
incardon committed
155 156
	Ghost<2,float> g(0.01);

157
	// Grid points on x=256 and y=64
158
	long int sz[] = {256,64};
159 160 161
	size_t szu[2];
	szu[0] = (size_t)sz[0];
	szu[1] = (size_t)sz[1];
Pietro Incardona's avatar
Pietro Incardona committed
162

163 164 165 166
	// We need one more point on the left and down part of the domain
	// This is given by the boundary conditions that we impose, the
	// reason is mathematical in order to have a well defined system
	// and cannot be discussed here
167
	Padding<2> pd({1,1},{0,0});
incardon's avatar
incardon committed
168

169
	// Distributed grid that store the solution
170
	grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> g_dist(szu,domain,g);
incardon's avatar
incardon committed
171

172
	// It is the maximum extension of the stencil
173 174
	Ghost<2,long int> stencil_max(1);

175
	// Finite difference scheme
176
	FDScheme<lid_nn> fd(pd, stencil_max, domain,g_dist);
177

178 179 180 181 182
	// Here we impose the equation, we start from the incompressibility Eq imposed in the bulk with the
	// exception of the first point {0,0} and than we set P = 0 in {0,0}, why we are doing this is again
	// mathematical to have a well defined system, an intuitive explanation is that P and P + c are both
	// solution for the incompressibility equation, this produce an ill-posed problem to make it well posed
	// we set one point in this case {0,0} the pressure to a fixed constant for convenience P = 0
183 184
	fd.impose(ic_eq(),0.0, EQ_3, {0,0},{sz[0]-2,sz[1]-2},true);
	fd.impose(Prs(),  0.0, EQ_3, {0,0},{0,0});
185 186

	// Here we impose the Eq1 and Eq2
187 188
	fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2});
	fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2});
189

190 191
	// v_x and v_y
	// Imposing B1
192 193
	fd.impose(v_x(),0.0, EQ_1, {0,0},{0,sz[1]-2});
	fd.impose(avg_vy_f(),0.0, EQ_2 , {-1,0},{-1,sz[1]-1});
194 195
	// Imposing B2
	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 199
	// Imposing B3
	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 202
	// Imposing B4
	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 206 207 208 209 210
	// When we pad the grid, there are points of the grid that are not
	// touched by the previous condition. Mathematically this lead
	// to have too many variables for the conditions that we are imposing.
	// Here we are imposing variables that we do not touch to zero
	//

211
	// Padding pressure
212 213 214 215
	fd.impose(Prs(), 0.0, EQ_3, {-1,-1},{sz[0]-1,-1});
	fd.impose(Prs(), 0.0, EQ_3, {-1,sz[1]-1},{sz[0]-1,sz[1]-1});
	fd.impose(Prs(), 0.0, EQ_3, {-1,0},{-1,sz[1]-2});
	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 219
	fd.impose(v_x(), 0.0, EQ_1, {-1,-1},{-1,sz[1]-1});
	fd.impose(v_y(), 0.0, EQ_2, {-1,-1},{sz[0]-1,-1});
220

Pietro Incardona's avatar
Pietro Incardona committed
221 222
	solver_type solver;
	auto x = solver.solve(fd.getA(),fd.getB());
223

224 225
	//! [lid-driven cavity 2D]

226 227
	//! [Copy the solution to grid]

Pietro Incardona's avatar
Pietro Incardona committed
228 229 230 231
	fd.template copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1},g_dist);

	std::string s = std::string(demangle(typeid(solver_type).name()));
	s += "_";
232 233 234

	//! [Copy the solution to grid]

Pietro Incardona's avatar
Pietro Incardona committed
235
	g_dist.write(s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid");
Pietro Incardona's avatar
Pietro Incardona committed
236

Pietro Incardona's avatar
Pietro Incardona committed
237 238
#ifdef HAVE_OSX

Pietro Incardona's avatar
Pietro Incardona committed
239 240
    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";
    std::string file2 = s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
Pietro Incardona's avatar
Pietro Incardona committed
241 242 243

#else

incardon's avatar
incardon committed
244
	#if __GNUC__ == 6
Pietro Incardona's avatar
Pietro Incardona committed
245

incardon's avatar
incardon committed
246 247 248 249 250 251 252
	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";
	std::string file2 = s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";

	#elif __GNUC__ == 5

	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";
	std::string file2 = s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
Pietro Incardona's avatar
Pietro Incardona committed
253 254 255

	#else

incardon's avatar
incardon committed
256 257
	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";
	std::string file2 = s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
Pietro Incardona's avatar
Pietro Incardona committed
258 259 260

	#endif

Pietro Incardona's avatar
Pietro Incardona committed
261

Pietro Incardona's avatar
Pietro Incardona committed
262 263
#endif

Pietro Incardona's avatar
Pietro Incardona committed
264 265 266
    std::cout << "File1: " << file1 << std::endl;
    std::cout << "File2: " << file2 << std::endl;

incardon's avatar
incardon committed
267 268
#ifndef SE_CLASS3

Pietro Incardona's avatar
Pietro Incardona committed
269 270 271
	// Check that match
	bool test = compare(file1,file2);
	BOOST_REQUIRE_EQUAL(test,true);
272

incardon's avatar
incardon committed
273 274
#endif

Pietro Incardona's avatar
Pietro Incardona committed
275 276 277 278 279 280 281
}

// Lid driven cavity, incompressible fluid

BOOST_AUTO_TEST_CASE(lid_driven_cavity)
{
	lid_driven_cavity_2d<umfpack_solver<double>,lid_nn>();
incardon's avatar
incardon committed
282
}
incardon's avatar
incardon committed
283

incardon's avatar
incardon committed
284 285
BOOST_AUTO_TEST_SUITE_END()

incardon's avatar
incardon committed
286
#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_HPP_ */