FDScheme.hpp 9.99 KB
Newer Older
incardon's avatar
incardon committed
1 2 3 4 5 6 7 8 9 10
/*
 * FiniteDifferences.hpp
 *
 *  Created on: Sep 17, 2015
 *      Author: i-bird
 */

#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_

11
#include "../Matrix/SparseMatrix.hpp"
incardon's avatar
incardon committed
12
#include "Grid/grid_dist_id.hpp"
incardon's avatar
incardon committed
13
#include "Grid/grid_dist_id_iterator_sub.hpp"
incardon's avatar
incardon committed
14
#include "eq.hpp"
15
#include "data_type/scalar.hpp"
16
#include "NN/CellList/CellDecomposer.hpp"
incardon's avatar
incardon committed
17 18

/*! \brief Finite Differences
19 20
 *
 * This class is able to discreatize on a Matrix any operator. In order to create a consistent
incardon's avatar
incardon committed
21 22
 * Matrix it is required that each processor must contain a contiguos range on grid points without
 * holes. In order to ensure this, each processor produce a contiguos local labelling of its local
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
 * points. Each processor also add an offset equal to the number of local
 * points of the processors with id smaller than him, to produce a global and non overlapping
 * labelling. An example is shown in the figures down, here we have
 * a grid 8x6 divided across two processor each processor label locally its grid points
 *
 * \verbatim
 *
+--------------------------+
| 1   2   3   4| 1  2  3  4|
|              |           |
| 5   6   7   8| 5  6  7  8|
|              |           |
| 9  10  11  12| 9 10 11 12|
+--------------------------+
|13  14  15| 13 14 15 16 17|
|          |               |
|16  17  18| 18 19 20 21 22|
|          |               |
|19  20  21| 23 24 25 26 27|
+--------------------------+

 *
 *
 * \endverbatim
 *
 * To the local relabelling is added an offset to make the local id global and non overlapping
 *
 *
 * \verbatim
 *
+--------------------------+
| 1   2   3   4|23 24 25 26|
|              |           |
| 5   6   7   8|27 28 29 30|
|              |           |
| 9  10  12  13|31 32 33 34|
+--------------------------+
|14  15  16| 35 36 37 38 39|
|          |               |
|17  18  19| 40 41 42 43 44|
|          |               |
|20  21  22| 45 46 47 48 49|
+--------------------------+
 *
 *
 * \endverbatim
incardon's avatar
incardon committed
69 70 71 72 73 74 75 76
 *
 * \tparam dim Dimensionality of the finite differences scheme
 *
 */

template<typename Sys_eqs>
class FDScheme
{
77 78 79 80 81 82 83 84 85
public:

	// Distributed grid map
	typedef grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition> g_map_type;

	typedef Sys_eqs Sys_eqs_typ;

private:

86 87 88
	// Padding
	Padding<Sys_eqs::dims> pd;

89 90
	typedef typename Sys_eqs::SparseMatrix_type::triplet_type triplet;

91
	// Sparse Matrix
92
	openfpm::vector<triplet> trpl;
incardon's avatar
incardon committed
93

94 95 96 97 98
	openfpm::vector<typename Sys_eqs::stype> b;

	// Domain Grid informations
	const grid_sm<Sys_eqs::dims,void> & gs;

99 100
	// Get the grid spacing
	typename Sys_eqs::stype spacing[Sys_eqs::dims];
101

102
	// mapping grid
103
	g_map_type g_map;
104 105 106 107 108 109 110

	// row of the matrix
	size_t row;

	// row on b
	size_t row_b;

111 112 113 114 115 116 117 118 119
	// Grid points that has each processor
	openfpm::vector<size_t> pnt;

	// Each point in the grid has a global id, to decompose correctly the Matrix each processor contain a
	// contiguos range of global id, example processor 0 can have from 0 to 234 and processor 1 from 235 to 512
	// no processors can have holes in the sequence, this number indicate where the sequence start for this
	// processor
	size_t s_pnt;

120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
	/* \brief calculate the mapping grid size with padding
	 *
	 * \param gs original grid size
	 *
	 * \return padded grid size
	 *
	 */
	inline const std::vector<size_t> padding( const size_t (& sz)[Sys_eqs::dims], Padding<Sys_eqs::dims> & pd)
	{
		std::vector<size_t> g_sz_pad(Sys_eqs::dims);

		for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
			g_sz_pad[i] = sz[i] + pd.getLow(i) + pd.getHigh(i);

		return g_sz_pad;
	}

137 138 139 140 141 142 143 144 145 146
	/*! \brief Check if the Matrix is consistent
	 *
	 */
	void consistency()
	{
		// A and B must have the same rows
		if (row != row_b)
			std::cerr << "Error " << __FILE__ << ":" << __LINE__ << "the term B and the Matrix A for Ax=B must contain the same number of rows\n";

		// Indicate all the non zero rows
147
		openfpm::vector<unsigned char> nz_rows;
148 149 150 151
		nz_rows.resize(row_b);

		for (size_t i = 0 ; i < trpl.size() ; i++)
		{
152
			nz_rows.get(trpl.get(i).row()) = true;
153 154 155
		}

		// Indicate all the non zero colums
156 157 158
		openfpm::vector<unsigned> nz_cols;
		nz_cols.resize(row_b);

159 160
		for (size_t i = 0 ; i < trpl.size() ; i++)
		{
161
			nz_cols.get(trpl.get(i).col()) = true;
162 163 164 165 166 167
		}

		// all the rows must have a non zero element
		for (size_t i = 0 ; i < nz_rows.size() ; i++)
		{
			if (nz_rows.get(i) == false)
168
				std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix row " << i <<  " is not filled\n";
169 170 171 172 173 174
		}

		// all the colums must have a non zero element
		for (size_t i = 0 ; i < nz_cols.size() ; i++)
		{
			if (nz_cols.get(i) == false)
175
				std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix colum " << i << " is not filled\n";
176 177 178
		}
	}

179
	/* \brief Convert discrete ghost into continous ghost
incardon's avatar
incardon committed
180
	 *
181
	 * \param Ghost
incardon's avatar
incardon committed
182 183
	 *
	 */
184
	inline Ghost<Sys_eqs::dims,typename Sys_eqs::stype> convert_dg_cg(const Ghost<Sys_eqs::dims,typename Sys_eqs::stype> & g)
incardon's avatar
incardon committed
185
	{
186
		return g_map_type::convert_ghost(g,g_map.getCellDecomposer());
incardon's avatar
incardon committed
187 188
	}

incardon's avatar
incardon committed
189 190
public:

191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
	/*! \brief Get the grid padding
	 *
	 * \return the grid padding
	 *
	 */
	const Padding<Sys_eqs::dims> & getPadding()
	{
		return pd;
	}

	/*! \brief Return the map between the grid index position and the position in the distributed vector
	 *
	 * \return the map
	 *
	 */
	const g_map_type & getMap()
	{
		return g_map;
	}

211 212 213 214 215 216
	/*! \brief Constructor
	 *
	 * \param pd Padding
	 * \param gs grid infos where Finite differences work
	 *
	 */
217 218
	FDScheme(Padding<Sys_eqs::dims> & pd, const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain, const grid_sm<Sys_eqs::dims,void> & gs, const typename Sys_eqs::b_grid::decomposition & dec, Vcluster & v_cl)
	:pd(pd),gs(gs),g_map(dec,padding(gs.getSize(),pd),domain,Ghost<Sys_eqs::dims,typename Sys_eqs::stype>(1)),row(0),row_b(0)
219
	{
220 221 222 223 224
		// Calculate the size of the local domain
		size_t sz = g_map.getLocalDomainSize();

		// Get the total size of the local grids on each processors
		v_cl.allGather(sz,pnt);
225 226
		v_cl.execute();
		s_pnt = 0;
227 228 229

		// calculate the starting point for this processor
		for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
230
			s_pnt += pnt.get(i);
231 232 233

		// Calculate the starting point

234 235 236 237 238 239 240 241 242 243
		// Counter
		size_t cnt = 0;

		// Create the re-mapping-grid
		auto it = g_map.getDomainIterator();

		while (it.isNext())
		{
			auto key = it.get();

244
			g_map.template get<0>(key) = cnt + s_pnt;
245 246 247 248 249

			++cnt;
			++it;
		}

250
		// sync the ghost
251

252
		g_map.template ghost_get<0>();
253 254 255 256 257 258 259 260 261 262 263

		// Create a CellDecomposer and calculate the spacing

		size_t sz_g[Sys_eqs::dims];
		for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
			sz_g[i] = gs.getSize()[i] - 1;

		CellDecomposer_sm<Sys_eqs::dims,typename Sys_eqs::stype> cd(domain,sz_g,0);

		for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
			spacing[i] = cd.getCellBox().getHigh(i);
264 265
	}

incardon's avatar
incardon committed
266 267
	/*! \brief Impose an operator
	 *
268 269 270 271 272
	 * This function impose an operator on a box region to produce the system
	 *
	 * Ax = b
	 *
	 * ## Stokes equation, lid driven cavity with one splipping wall
incardon's avatar
incardon committed
273
	 *
274 275 276 277 278
	 * \param op Operator to impose (A term)
	 * \param num right hand side of the term (b term)
	 * \param id Equation id in the system that we are imposing
	 * \param start starting point of the box
	 * \param stop stop point of the box
incardon's avatar
incardon committed
279 280
	 *
	 */
281
	template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,const long int (& start)[Sys_eqs::dims], const long int (& stop)[Sys_eqs::dims], bool skip_first = false)
incardon's avatar
incardon committed
282
	{
283
		// add padding to start and stop
284 285
		grid_key_dx<Sys_eqs::dims> start_k = grid_key_dx<Sys_eqs::dims>(start) + pd.getKP1();
		grid_key_dx<Sys_eqs::dims> stop_k = grid_key_dx<Sys_eqs::dims>(stop) + pd.getKP1();
286

287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309
		auto it = g_map.getSubDomainIterator(start_k,stop_k);

		impose(op,num,id,it,skip_first);
	}

	/*! \brief Impose an operator
	 *
	 * This function impose an operator on a particular grid region to produce the system
	 *
	 * Ax = b
	 *
	 * ## Stokes equation, lid driven cavity with one splipping wall
	 *
	 * \param op Operator to impose (A term)
	 * \param num right hand side of the term (b term)
	 * \param id Equation id in the system that we are imposing
	 * \param it_d iterator that define where you want to impose
	 *
	 */
	template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d, bool skip_first = false)
	{
		auto it = it_d;
		grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
310

incardon's avatar
incardon committed
311 312
		std::unordered_map<long int,float> cols;

313 314 315
		// resize b if needed
		b.resize(Sys_eqs::nvar * gs.size());

316 317
		bool is_first = skip_first;

incardon's avatar
incardon committed
318 319 320
		// iterate all the grid points
		while (it.isNext())
		{
321 322 323 324 325 326 327
			if (is_first == true)
			{
				++it;
				is_first = false;
				continue;
			}

incardon's avatar
incardon committed
328 329
			// get the position
			auto key = it.get();
330
			grid_key_dx<2> gkey = g_map.getGKey(key);
incardon's avatar
incardon committed
331

incardon's avatar
incardon committed
332
			// Calculate the non-zero colums
333
			T::value(g_map,key,gs,spacing,cols,1.0);
incardon's avatar
incardon committed
334 335 336 337 338 339

			// create the triplet

			for ( auto it = cols.begin(); it != cols.end(); ++it )
			{
				trpl.add();
340
				trpl.last().row() = Sys_eqs::nvar * gs.LinId(gkey) + id;
341 342
				trpl.last().col() = it->first;
				trpl.last().value() = it->second;
incardon's avatar
incardon committed
343

344
				std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
345
			}
incardon's avatar
incardon committed
346

347 348
			b.get(Sys_eqs::nvar * gs.LinId(gkey) + id) = num;

incardon's avatar
incardon committed
349
			cols.clear();
350
			std::cout << "\n";
incardon's avatar
incardon committed
351

352 353 354 355
			++row;
			++row_b;
			++it;
		}
incardon's avatar
incardon committed
356
	}
incardon's avatar
incardon committed
357

358 359
	typename Sys_eqs::SparseMatrix_type A;

incardon's avatar
incardon committed
360 361 362 363 364
	/*! \brief produce the Matrix
	 *
	 *  \tparam Syst_eq System of equations, or a single equation
	 *
	 */
365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384
	typename Sys_eqs::SparseMatrix_type & getA()
	{
#ifdef SE_CLASS1
		consistency();
#endif
		A.resize(row,row);
		A.setFromTriplets(trpl);

		return A;

	}

	typename Sys_eqs::Vector_type B;

	/*! \brief produce the B vector
	 *
	 *  \tparam Syst_eq System of equations, or a single equation
	 *
	 */
	typename Sys_eqs::Vector_type & getB()
incardon's avatar
incardon committed
385
	{
386 387 388
#ifdef SE_CLASS1
		consistency();
#endif
incardon's avatar
incardon committed
389

390 391 392 393 394
		B.resize(row_b);

		// copy the vector
		for (size_t i = 0; i < row_b; i++)
			B.get(i) = b.get(i);
incardon's avatar
incardon committed
395

396
		return B;
incardon's avatar
incardon committed
397 398 399 400 401 402 403 404
	}
};

#define EQS_FIELDS 0
#define EQS_SPACE 1


#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_ */