FDScheme.hpp 15.8 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"
17
#include "Grid/staggered_dist_grid_util.hpp"
incardon's avatar
incardon committed
18 19

/*! \brief Finite Differences
20
 *
Pietro Incardona's avatar
Pietro Incardona committed
21 22 23
 * This class is able to discretize on a Matrix any system of equations producing a linear system of type \f$Ax=b\f$. In order to create a consistent
 * Matrix it is required that each processor must contain a contiguous range on grid points without
 * holes. In order to ensure this, each processor produce a contiguous local labeling of its local
24 25
 * 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
Pietro Incardona's avatar
Pietro Incardona committed
26 27
 * labeling. An example is shown in the figures down, here we have
 * a grid 8x6 divided across four processors each processor label locally its grid points
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 69
 *
 * \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
70
 *
71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
 * \tparam Sys_eqs Definition of the system of equations
 *
 * # Examples
 *
 * ## Solve lid-driven cavity 2D for incompressible fluid (inertia=0 --> Re=0)
 *
 * In this case the system of equation to solve is
 *
 * \f$
\left\{
\begin{array}{c}
\eta\nabla v_x + \partial_x P = 0 \quad Eq1 \\
\eta\nabla v_y + \partial_y P = 0 \quad Eq2 \\
\partial_x v_x + \partial_y v_y = 0 \quad Eq3
\end{array}
\right.  \f$

and  boundary conditions

 * \f$
\left\{
\begin{array}{c}
v_x = 0, v_y = 0 \quad x = 0 \quad B1\\
v_x = 0, v_y = 1.0 \quad x = L \quad B2\\
v_x = 0, v_y = 0 \quad y = 0 \quad B3\\
v_x = 0, v_y = 0 \quad y = L \quad B4\\
\end{array}
\right.  \f$

 *
 * with \f$v_x\f$ and \f$v_y\f$ the velocity in x and y and \f$P\f$ Pressure
 *
 * In order to solve such system first we define the general properties of the system
 *
 *	\snippet eq_unit_test.hpp Definition of the system
 *
 * ## Define the equations of the system
 *
 * \snippet eq_unit_test.hpp Definition of the equation of the system in the bulk and at the boundary
 *
 * ## Define the domain and impose the equations
 *
 * \snippet eq_unit_test.hpp lid-driven cavity 2D
 *
 * # 3D
 *
 * A 3D case is given in the examples
incardon's avatar
incardon committed
118 119 120 121 122 123
 *
 */

template<typename Sys_eqs>
class FDScheme
{
124 125 126
public:

	// Distributed grid map
Pietro Incardona's avatar
Pietro Incardona committed
127
	typedef grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition::extended_type> g_map_type;
128 129 130 131 132

	typedef Sys_eqs Sys_eqs_typ;

private:

133 134 135
	// Padding
	Padding<Sys_eqs::dims> pd;

136 137
	typedef typename Sys_eqs::SparseMatrix_type::triplet_type triplet;

Pietro Incardona's avatar
Pietro Incardona committed
138 139
	// Vector b
	typename Sys_eqs::Vector_type b;
140 141 142 143

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

144 145
	// Get the grid spacing
	typename Sys_eqs::stype spacing[Sys_eqs::dims];
146

147
	// mapping grid
148
	g_map_type g_map;
149 150 151 152 153 154 155

	// row of the matrix
	size_t row;

	// row on b
	size_t row_b;

156 157 158
	// Grid points that has each processor
	openfpm::vector<size_t> pnt;

159 160 161
	// Staggered position for each property
	comb<Sys_eqs::dims> s_pos[Sys_eqs::nvar];

162 163 164 165 166 167
	// 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;

168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
	struct key_and_eq
	{
		grid_key_dx<Sys_eqs::dims> key;
		size_t eq;
	};

	/*! \brief From the row Matrix position to the spatial position
	 *
	 * \param row Matrix
	 *
	 * \return spatial position
	 *
	 */
	inline key_and_eq from_row_to_key(size_t row)
	{
		key_and_eq ke;
		auto it = g_map.getDomainIterator();

		while (it.isNext())
		{
			size_t row_low = g_map.template get<0>(it.get());

			if (row >= row_low * Sys_eqs::nvar && row < row_low * Sys_eqs::nvar + Sys_eqs::nvar)
			{
				ke.eq = row - row_low * Sys_eqs::nvar;
				ke.key = g_map.getGKey(it.get());
				ke.key -= pd.getKP1();
				return ke;
			}

			++it;
		}
		std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " the row does not map to any position" << "\n";

		return ke;
	}

205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
	/* \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;
	}

222 223 224 225 226
	/*! \brief Check if the Matrix is consistent
	 *
	 */
	void consistency()
	{
227 228
		openfpm::vector<triplet> & trpl = A.getMatrixTriplets();

229 230 231 232 233
		// 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
234
		openfpm::vector<unsigned char> nz_rows;
235 236 237 238
		nz_rows.resize(row_b);

		for (size_t i = 0 ; i < trpl.size() ; i++)
		{
239
			nz_rows.get(trpl.get(i).row()) = true;
240 241 242
		}

		// Indicate all the non zero colums
243 244 245
		openfpm::vector<unsigned> nz_cols;
		nz_cols.resize(row_b);

246 247
		for (size_t i = 0 ; i < trpl.size() ; i++)
		{
248
			nz_cols.get(trpl.get(i).col()) = true;
249 250 251 252 253 254
		}

		// 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)
255 256 257 258
			{
				key_and_eq ke = from_row_to_key(i);
				std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix row " << i <<  " is not filled, position " << ke.key.to_string() << " equation: " << ke.eq << "\n";
			}
259 260 261 262 263 264
		}

		// 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)
265
				std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix colum " << i << " is not filled\n";
266 267 268
		}
	}

Pietro Incardona's avatar
Pietro Incardona committed
269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318
	/*! \brief Copy a given solution vector in a staggered grid
	 *
	 * \tparam Vct Vector type
	 * \tparam Grid_dst target grid
	 * \tparam pos set of properties
	 *
	 * \param v Vector
	 * \param g_dst target staggered grid
	 *
	 */
	template<typename Vct, typename Grid_dst ,unsigned int ... pos> void copy_staggered(Vct & v, Grid_dst & g_dst)
	{
		// check that g_dst is staggered
		if (g_dst.is_staggered() == false)
			std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid must be staggered " << std::endl;

#ifdef SE_CLASS1

		if (g_map.getLocalDomainSize() != g_dst.getLocalDomainSize())
			std::cerr << __FILE__ << ":" << __LINE__ << " The staggered and destination grid in size does not match " << std::endl;
#endif

		// sub-grid iterator over the grid map
		auto g_map_it = g_map.getDomainIterator();

		// Iterator over the destination grid
		auto g_dst_it = g_dst.getDomainIterator();

		while (g_map_it.isNext() == true)
		{
			typedef typename to_boost_vmpl<pos...>::type vid;
			typedef boost::mpl::size<vid> v_size;

			auto key_src = g_map_it.get();
			size_t lin_id = g_map.template get<0>(key_src);

			// destination point
			auto key_dst = g_dst_it.get();

			// Transform this id into an id for the Eigen vector

			copy_ele<Sys_eqs_typ, Grid_dst,Vct> cp(key_dst,g_dst,v,lin_id,g_map.size());

			boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);

			++g_map_it;
			++g_dst_it;
		}
	}

319 320 321
public:

	/*! \brief set the staggered position for each property
incardon's avatar
incardon committed
322
	 *
323
	 * \param vector containing the staggered position for each property
incardon's avatar
incardon committed
324 325
	 *
	 */
326
	void setStagPos(comb<Sys_eqs::dims> (& sp)[Sys_eqs::nvar])
incardon's avatar
incardon committed
327
	{
328 329
		for (size_t i = 0 ; i < Sys_eqs::nvar ; i++)
			s_pos[i] = sp[i];
incardon's avatar
incardon committed
330 331
	}

332 333 334 335 336 337 338 339 340 341 342 343
	/*! \brief compute the staggered position for each property
	 *
	 * This is compute from the value_type stored by Sys_eqs::b_grid::value_type
	 *
	 * ### Example
	 *
	 * \snippet eq_unit_test.hpp Compute staggered properties
	 *
	 */
	void computeStag()
	{
		typedef typename Sys_eqs::b_grid::value_type prp_type;
incardon's avatar
incardon committed
344

345 346 347 348 349 350 351 352
		openfpm::vector<comb<Sys_eqs::dims>> c_prp[prp_type::max_prop];

		stag_set_position<Sys_eqs::dims,prp_type> ssp(c_prp);

		boost::mpl::for_each_ref< boost::mpl::range_c<int,0,prp_type::max_prop> >(ssp);
	}

	/*! \brief Get the specified padding
353
	 *
354
	 * \return the padding specified
355 356 357 358 359 360 361 362
	 *
	 */
	const Padding<Sys_eqs::dims> & getPadding()
	{
		return pd;
	}

	/*! \brief Return the map between the grid index position and the position in the distributed vector
363 364
	 *
	 * It is the map explained in the intro of the FDScheme
365 366 367 368 369 370 371 372 373
	 *
	 * \return the map
	 *
	 */
	const g_map_type & getMap()
	{
		return g_map;
	}

374 375
	/*! \brief Constructor
	 *
376 377
	 * \param pd Padding, how many points out of boundary are present
	 * \param domain extension of the domain
378
	 * \param gs grid infos where Finite differences work
Pietro Incardona's avatar
Pietro Incardona committed
379
	 * \param stencil maximum extension of the stencil on each directions
380
	 * \param dec Decomposition of the domain
381 382
	 *
	 */
Pietro Incardona's avatar
Pietro Incardona committed
383 384
	FDScheme(Padding<Sys_eqs::dims> & pd, const Ghost<Sys_eqs::dims,long int> & stencil, const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain, const grid_sm<Sys_eqs::dims,void> & gs, const typename Sys_eqs::b_grid & b_g)
	:pd(pd),gs(gs),g_map(b_g,stencil,pd),row(0),row_b(0)
385
	{
Pietro Incardona's avatar
Pietro Incardona committed
386
		Vcluster & v_cl = b_g.getDecomposition().getVC();
387

388 389 390 391 392
		// 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);
393 394
		v_cl.execute();
		s_pnt = 0;
395 396 397

		// calculate the starting point for this processor
		for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
398
			s_pnt += pnt.get(i);
399 400 401

		// Calculate the starting point

402 403 404
		// Counter
		size_t cnt = 0;

405
		// Create the re-mapping grid
406 407 408 409 410 411
		auto it = g_map.getDomainIterator();

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

412
			g_map.template get<0>(key) = cnt + s_pnt;
413 414 415 416 417

			++cnt;
			++it;
		}

418 419
		// sync the ghost
		g_map.template ghost_get<0>();
420 421 422 423 424 425 426 427 428 429 430

		// 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);
431 432
	}

incardon's avatar
incardon committed
433 434
	/*! \brief Impose an operator
	 *
435 436 437 438 439
	 * This function impose an operator on a box region to produce the system
	 *
	 * Ax = b
	 *
	 * ## Stokes equation, lid driven cavity with one splipping wall
Pietro Incardona's avatar
Pietro Incardona committed
440
	 * \snippet eq_unit_test.hpp lid-driven cavity 2D
incardon's avatar
incardon committed
441
	 *
442 443 444 445 446
	 * \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
447 448
	 *
	 */
449
	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
450
	{
451
		// add padding to start and stop
Pietro Incardona's avatar
Pietro Incardona committed
452 453
		grid_key_dx<Sys_eqs::dims> start_k = grid_key_dx<Sys_eqs::dims>(start);
		grid_key_dx<Sys_eqs::dims> stop_k = grid_key_dx<Sys_eqs::dims>(stop);
454

455 456 457 458 459 460 461 462 463 464 465
		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
	 *
Pietro Incardona's avatar
Pietro Incardona committed
466 467
	 * ## Stokes equation 2D, lid driven cavity with one splipping wall
	 * \snippet eq_unit_test.hpp Copy the solution to grid
468 469 470 471 472 473 474 475 476
	 *
	 * \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)
	{
477
		Vcluster & v_cl = create_vcluster();
Pietro Incardona's avatar
Pietro Incardona committed
478

479 480
		openfpm::vector<triplet> & trpl = A.getMatrixTriplets();

481 482
		auto it = it_d;
		grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
483

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

486
		// resize b if needed
Pietro Incardona's avatar
Pietro Incardona committed
487
		b.resize(Sys_eqs::nvar * g_map.size());
488

489 490
		bool is_first = skip_first;

incardon's avatar
incardon committed
491 492 493
		// iterate all the grid points
		while (it.isNext())
		{
Pietro Incardona's avatar
Pietro Incardona committed
494
			if (is_first == true && v_cl.getProcessUnitID() == 0)
495 496 497 498 499
			{
				++it;
				is_first = false;
				continue;
			}
Pietro Incardona's avatar
Pietro Incardona committed
500 501
			else
				is_first = false;
502

incardon's avatar
incardon committed
503 504 505
			// get the position
			auto key = it.get();

incardon's avatar
incardon committed
506
			// Calculate the non-zero colums
507
			T::value(g_map,key,gs,spacing,cols,1.0);
incardon's avatar
incardon committed
508 509 510 511 512 513

			// create the triplet

			for ( auto it = cols.begin(); it != cols.end(); ++it )
			{
				trpl.add();
514
				trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
515 516
				trpl.last().col() = it->first;
				trpl.last().value() = it->second;
incardon's avatar
incardon committed
517

518
//				std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
519
			}
incardon's avatar
incardon committed
520

Pietro Incardona's avatar
Pietro Incardona committed
521
			b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num;
522

incardon's avatar
incardon committed
523
			cols.clear();
524 525 526 527
//			std::cout << "\n";

			// if SE_CLASS1 is defined check the position
#ifdef SE_CLASS1
528
//			T::position(key,gs,s_pos);
529
#endif
incardon's avatar
incardon committed
530

531 532 533 534
			++row;
			++row_b;
			++it;
		}
incardon's avatar
incardon committed
535
	}
incardon's avatar
incardon committed
536

537 538
	typename Sys_eqs::SparseMatrix_type A;

incardon's avatar
incardon committed
539 540
	/*! \brief produce the Matrix
	 *
541
	 *  \return the Sparse matrix produced
incardon's avatar
incardon committed
542 543
	 *
	 */
544 545 546 547 548
	typename Sys_eqs::SparseMatrix_type & getA()
	{
#ifdef SE_CLASS1
		consistency();
#endif
Pietro Incardona's avatar
Pietro Incardona committed
549
		A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*Sys_eqs::nvar);
550 551 552 553 554 555 556

		return A;

	}

	/*! \brief produce the B vector
	 *
557
	 *  \return the vector produced
558 559 560
	 *
	 */
	typename Sys_eqs::Vector_type & getB()
incardon's avatar
incardon committed
561
	{
562 563 564
#ifdef SE_CLASS1
		consistency();
#endif
incardon's avatar
incardon committed
565

Pietro Incardona's avatar
Pietro Incardona committed
566
		return b;
incardon's avatar
incardon committed
567
	}
Pietro Incardona's avatar
Pietro Incardona committed
568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603

	/*! \brief Copy the vector into the grid
	 *
	 * ## Copy the solution into the grid
	 * \snippet eq_unit_test.hpp Copy the solution to grid
	 *
	 * \tparam scheme Discretization scheme
	 * \tparam Grid_dst type of the target grid
	 * \tparam pos target properties
	 *
	 * \param scheme Discretization scheme
	 * \param start point
	 * \param stop point
	 * \param g_dst Destination grid
	 *
	 */
	template<unsigned int ... pos, typename Vct, typename Grid_dst> void copy(Vct & v,const long int (& start)[Sys_eqs_typ::dims], const long int (& stop)[Sys_eqs_typ::dims], Grid_dst & g_dst)
	{
		if (is_grid_staggered<Sys_eqs>::value())
		{
			if (g_dst.is_staggered() == true)
				copy_staggered<Vct,Grid_dst,pos...>(v,g_dst);
			else
			{
				// Create a temporal staggered grid and copy the data there
				auto & g_map = this->getMap();
				staggered_grid_dist<Grid_dst::dims,typename Grid_dst::stype,typename Grid_dst::value_type,typename Grid_dst::decomposition::extended_type, typename Grid_dst::memory_type, typename Grid_dst::device_grid_type> stg(g_dst,g_map.getDecomposition().getGhost(),this->getPadding());
				stg.setDefaultStagPosition();
				copy_staggered<Vct,decltype(stg),pos...>(v,stg);

				// sync the ghost and interpolate to the normal grid
				stg.template ghost_get<pos...>();
				stg.template to_normal<Grid_dst,pos...>(g_dst,this->getPadding(),start,stop);
			}
		}
	}
incardon's avatar
incardon committed
604 605 606 607 608 609 610
};

#define EQS_FIELDS 0
#define EQS_SPACE 1


#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_ */