eq_unit_test.hpp 7.95 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 17
#include "Grid/grid_dist_id.hpp"
#include "data_type/scalar.hpp"
#include "Decomposition/CartDecomposition.hpp"
18 19 20
#include "Vector/Vector.hpp"
#include "Solvers/umfpack_solver.hpp"
#include "data_type/aggregate.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
21
#include "FiniteDifference/FDScheme.hpp"
incardon's avatar
incardon committed
22

23 24
BOOST_AUTO_TEST_SUITE( eq_test_suite )

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

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

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

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

	// type of space float, double, ...
	typedef float stype;
incardon's avatar
incardon committed
40

41
	// type of base grid, it is the distributed grid that will store the result
Pietro Incardona's avatar
Pietro Incardona committed
42
	// Note the first property is a 2D vector (velocity), the second is a scalar (Pressure)
43 44
	typedef grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> b_grid;

45
	// type of SparseMatrix, for the linear system, this parameter is bounded by the solver
Pietro Incardona's avatar
Pietro Incardona committed
46
	// that you are using, in case of umfpack it is the only possible choice
47 48
	typedef SparseMatrix<double,int> SparseMatrix_type;

49
	// type of Vector for the linear system, this parameter is bounded by the solver
Pietro Incardona's avatar
Pietro Incardona committed
50
	// that you are using, in case of umfpack it is the only possible choice
51 52
	typedef Vector<double> Vector_type;

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

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

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

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

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

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

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

76
// Create field that we have v_x, v_y, P
incardon's avatar
incardon committed
77 78 79 80 81 82 83 84
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;
85 86
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
87 88 89 90

// Eq2 V_y

typedef mul<eta,Lap<v_y,lid_nn>,lid_nn> eta_lap_vy;
91 92 93
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
94 95 96

// Eq3 Incompressibility

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

incardon's avatar
incardon committed
101

102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
// 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
 *
 */

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

130 131 132 133 134 135 136
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

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

// Lid driven cavity, incompressible fluid
incardon's avatar
incardon committed
140

141
BOOST_AUTO_TEST_CASE(lid_driven_cavity)
incardon's avatar
incardon committed
142
{
143
	Vcluster & v_cl = create_vcluster();
Pietro Incardona's avatar
Pietro Incardona committed
144

Pietro Incardona's avatar
Pietro Incardona committed
145 146 147
	if (v_cl.getProcessingUnits() > 3)
		return;

148 149
	//! [lid-driven cavity 2D]

Pietro Incardona's avatar
Pietro Incardona committed
150 151 152 153
	// velocity in the grid is the property 0, pressure is the property 1
	constexpr int velocity = 0;
	constexpr int pressure = 1;

154
	// Domain, a rectangle
Pietro Incardona's avatar
Pietro Incardona committed
155
	Box<2,float> domain({0.0,0.0},{3.0,1.0});
incardon's avatar
incardon committed
156

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

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

166 167 168 169
	// 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
170
	Padding<2> pd({1,1},{0,0});
incardon's avatar
incardon committed
171

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

Pietro Incardona's avatar
Pietro Incardona committed
175
	// It is the maximum extension of the stencil
Pietro Incardona's avatar
Pietro Incardona committed
176 177
	Ghost<2,long int> stencil_max(1);

178
	// Finite difference scheme
Pietro Incardona's avatar
Pietro Incardona committed
179
	FDScheme<lid_nn> fd(pd, stencil_max, domain, g_dist.getGridInfo(), g_dist);
180

181 182 183 184 185
	// 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
186 187
	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});
188 189

	// Here we impose the Eq1 and Eq2
190 191
	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});
Pietro Incardona's avatar
Pietro Incardona committed
192

193 194
	// v_x and v_y
	// Imposing B1
195 196
	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});
197 198
	// Imposing B2
	fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
199
	fd.impose(avg_vy(),1.0, EQ_2,    {sz[0]-1,0},{sz[0]-1,sz[1]-1});
Pietro Incardona's avatar
Pietro Incardona committed
200

201 202
	// Imposing B3
	fd.impose(avg_vx_f(),0.0, EQ_1, {0,-1},{sz[0]-1,-1});
203
	fd.impose(v_y(), 0.0, EQ_2, {0,0},{sz[0]-2,0});
204 205
	// Imposing B4
	fd.impose(avg_vx(),0.0, EQ_1,   {0,sz[1]-1},{sz[0]-1,sz[1]-1});
206
	fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1},{sz[0]-2,sz[1]-1});
Pietro Incardona's avatar
Pietro Incardona committed
207

208 209 210 211 212 213
	// 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
	//

Pietro Incardona's avatar
Pietro Incardona committed
214
	// Padding pressure
215 216 217 218
	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});
Pietro Incardona's avatar
Pietro Incardona committed
219 220

	// Impose v_x Padding Impose v_y padding
221 222
	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});
223 224 225

	auto x = umfpack_solver<double>::solve(fd.getA(),fd.getB());

226 227
	//! [lid-driven cavity 2D]

Pietro Incardona's avatar
Pietro Incardona committed
228 229 230 231 232 233
	//! [Copy the solution to grid]

	fd.copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1},g_dist);

	//! [Copy the solution to grid]

Pietro Incardona's avatar
Pietro Incardona committed
234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262
	g_dist.write("lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()));

	if (v_cl.getProcessUnitID() == 0)
	{
		if (v_cl.getProcessingUnits() == 1)
		{
			// Check that match
			bool test = compare("lid_driven_cavity_p1_grid_0_test.vtk","lid_driven_cavity_grid_0.vtk");
			BOOST_REQUIRE_EQUAL(test,true);
		}
		else if (v_cl.getProcessingUnits() == 2)
		{
			// Check that match
			bool test = compare("lid_driven_cavity_p2_grid_0_test.vtk","lid_driven_cavity_p2_grid_0.vtk");
			BOOST_REQUIRE_EQUAL(test,true);
			test = compare("lid_driven_cavity_p2_grid_1_test.vtk","lid_driven_cavity_p2_grid_1.vtk");
			BOOST_REQUIRE_EQUAL(test,true);
		}
		else if (v_cl.getProcessingUnits() == 3)
		{
			// Check that match
			bool test = compare("lid_driven_cavity_p3_grid_0_test.vtk","lid_driven_cavity_p3_grid_0.vtk");
			BOOST_REQUIRE_EQUAL(test,true);
			test = compare("lid_driven_cavity_p3_grid_1_test.vtk","lid_driven_cavity_p3_grid_1.vtk");
			BOOST_REQUIRE_EQUAL(test,true);
			test = compare("lid_driven_cavity_p3_grid_2_test.vtk","lid_driven_cavity_p3_grid_2.vtk");
			BOOST_REQUIRE_EQUAL(test,true);
		}
	}
incardon's avatar
incardon committed
263
}
incardon's avatar
incardon committed
264

incardon's avatar
incardon committed
265 266
BOOST_AUTO_TEST_SUITE_END()

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