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

#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_

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

/*! \brief Finite Differences
18 19
 *
 * This class is able to discreatize on a Matrix any operator. In order to create a consistent
incardon's avatar
incardon committed
20 21
 * 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
22 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
 * 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
68 69 70 71 72 73 74 75
 *
 * \tparam dim Dimensionality of the finite differences scheme
 *
 */

template<typename Sys_eqs>
class FDScheme
{
76 77 78 79
	// Padding
	Padding<Sys_eqs::dims> pd;

	// Sparse Matrix
incardon's avatar
incardon committed
80 81
	openfpm::vector<triplet<typename Sys_eqs::stype>> trpl;

82 83 84 85 86 87
	openfpm::vector<typename Sys_eqs::stype> b;

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

	// mapping grid
88
	grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition> g_map;
89 90 91 92 93 94 95

	// row of the matrix
	size_t row;

	// row on b
	size_t row_b;

96 97 98 99 100 101 102 103 104
	// 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;

105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
	/*! \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
		openfpm::vector<bool> nz_rows;
		nz_rows.resize(row_b);

		for (size_t i = 0 ; i < trpl.size() ; i++)
		{
			nz_rows.get(trpl.get(i).i) = true;
		}

		// Indicate all the non zero colums
		openfpm::vector<bool> nz_cols;
		for (size_t i = 0 ; i < trpl.size() ; i++)
		{
			nz_cols.get(trpl.get(i).j) = true;
		}

		// 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)
				std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix not all the row are filled\n";
		}

		// 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)
				std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix not all the row are filled\n";
		}
	}

incardon's avatar
incardon committed
145 146 147 148 149 150 151 152 153 154
	/*! \brief Convert from integer ghost to continuos
	 *
	 * \return the continuos version of the ghost
	 *
	 */
	Ghost<dim,Sys_eqs::stype> convert_into_cg()
	{

	}

incardon's avatar
incardon committed
155 156
public:

157 158 159 160 161 162
	/*! \brief Constructor
	 *
	 * \param pd Padding
	 * \param gs grid infos where Finite differences work
	 *
	 */
163
	FDScheme(Padding<Sys_eqs::dims> & pd, const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::b_grid::decomposition & dec, Vcluster & v_cl)
incardon's avatar
incardon committed
164
	:pd(pd),gs(gs),g_map(dec.duplicate(Ghost<Sys_eqs::dims,long int>(1)),gs.getSize(),domain)
165
	{
166 167 168 169 170 171 172 173 174 175 176 177
		// 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);

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

		// Calculate the starting point

178 179 180 181 182 183 184 185 186 187 188 189 190
		// Counter
		size_t cnt = 0;

		// Resize the b term
		b.resize(gs.size());

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

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

191
			g_map.template get<0>(key) = cnt + s_pnt;
192 193 194 195 196

			++cnt;
			++it;
		}

197
		// sync the ghost
198

199
		g_map.template ghost_get<0>();
200 201
	}

incardon's avatar
incardon committed
202 203 204 205 206
	/*! \brief Impose an operator
	 *
	 *
	 *
	 */
207
	template<typename T> void imposeA(const T & op , grid_dist_iterator_sub<Sys_eqs::dims,typename Sys_eqs::b_grid::d_grid> it)
incardon's avatar
incardon committed
208 209 210 211 212 213 214 215 216
	{
		std::unordered_map<long int,float> cols;

		// iterate all the grid points
		while (it.isNext())
		{
			// get the position
			auto key = it.get();

incardon's avatar
incardon committed
217
			// Calculate the non-zero colums
218
			T::value(g_map,key,gs,cols,1.0);
incardon's avatar
incardon committed
219 220 221 222 223 224

			// create the triplet

			for ( auto it = cols.begin(); it != cols.end(); ++it )
			{
				trpl.add();
225
				trpl.last().i = row;
incardon's avatar
incardon committed
226 227 228
				trpl.last().j = it->first;
				trpl.last().value = it->second;

229 230
				std::cout << "(" << trpl.last().i << "," << trpl.last().j << "," << trpl.last().value << ")" << "\n";
			}
incardon's avatar
incardon committed
231 232

			cols.clear();
233
			std::cout << "\n";
incardon's avatar
incardon committed
234

235
			++row;
incardon's avatar
incardon committed
236 237
			++it;
		}
incardon's avatar
incardon committed
238 239
	}

240 241
	/*! \brief Impose an operator
	 *
incardon's avatar
incardon committed
242 243 244
	 *
	 *
	 */
245
	void imposeB(typename Sys_eqs::stype num , grid_dist_iterator_sub<Sys_eqs::dims,typename Sys_eqs::b_grid::d_grid> it)
incardon's avatar
incardon committed
246
	{
247 248 249 250 251
		std::unordered_map<long int,float> cols;

		// iterate all the grid points
		while (it.isNext())
		{
incardon's avatar
incardon committed
252

253
			b.get(row_b) = num;
incardon's avatar
incardon committed
254

255 256 257
			++row_b;
			++it;
		}
incardon's avatar
incardon committed
258
	}
incardon's avatar
incardon committed
259 260 261 262 263 264

	/*! \brief produce the Matrix
	 *
	 *  \tparam Syst_eq System of equations, or a single equation
	 *
	 */
265
	template<typename Grid> SparseMatrix<typename Sys_eqs::stype> getMatrix()
incardon's avatar
incardon committed
266
	{
267 268 269
#ifdef SE_CLASS1
		consistency();
#endif
incardon's avatar
incardon committed
270 271 272 273 274 275 276 277 278 279


	}
};

#define EQS_FIELDS 0
#define EQS_SPACE 1


#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_ */