Derivative.hpp 13.8 KB
Newer Older
incardon's avatar
incardon committed
1 2 3 4
/*
 * Derivative.hpp
 *
 *  Created on: Oct 5, 2015
5
 *      Author: Pietro Incardona
incardon's avatar
incardon committed
6 7 8 9 10 11 12
 */

#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_

#define CENTRAL 0
#define CENTRAL_B_ONE_SIDE 1
13
#define FORWARD 2
14
#define BACKWARD 3
incardon's avatar
incardon committed
15 16 17 18 19

#include "util/mathutil.hpp"
#include "Vector/map_vector.hpp"
#include "Grid/comb.hpp"
#include "FiniteDifference/util/common.hpp"
incardon's avatar
incardon committed
20
#include "util/util_num.hpp"
21
#include <unordered_map>
22

incardon's avatar
incardon committed
23 24 25 26 27 28 29 30 31 32
/*! \brief Derivative second order on h (spacing)
 *
 * \tparam d on which dimension derive
 * \tparam Field which field derive
 * \tparam impl which implementation
 *
 */
template<unsigned int d, typename Field, typename Sys_eqs, unsigned int impl=CENTRAL>
class D
{
33 34 35 36 37 38
	/*! \brief Calculate which colums of the Matrix are non zero
	 *
	 * \param pos position where the derivative is calculated
	 * \param gs Grid info
	 * \param cols non-zero colums calculated by the function
	 * \param coeff coefficent (constant in front of the derivative)
incardon's avatar
incardon committed
39
	 *
40 41 42
	 * ### Example
	 *
	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
incardon's avatar
incardon committed
43 44
	 *
	 */
incardon's avatar
incardon committed
45
	inline static void value(const grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
incardon's avatar
incardon committed
46 47 48 49 50 51
	{
		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
	}

	/*! \brief Calculate the position where the derivative is calculated
	 *
52 53 54
	 * In case of non staggered case this function just return a null grid_key, in case of staggered,
	 *  it calculate how the operator shift the calculation in the cell
	 *
55
	 * \param pos position where we are calculating the derivative
56 57
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
incardon's avatar
incardon committed
58
	 *
59 60
	 * \return where (in which cell) the derivative is calculated
	 *
incardon's avatar
incardon committed
61
	 */
62 63 64
	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos,
			                                          const grid_sm<Sys_eqs::dims,void> & gs,
													  const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
incardon's avatar
incardon committed
65 66
	{
		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
67 68

		return pos;
incardon's avatar
incardon committed
69 70 71
	}
};

72
/*! \brief Second order central Derivative scheme on direction i
incardon's avatar
incardon committed
73
 *
74 75 76 77 78 79
 * \verbatim
 *
 *  -1        +1
 *   *---+---*
 *
 * \endverbatim
incardon's avatar
incardon committed
80 81 82 83 84 85 86
 *
 */
template<unsigned int d, typename arg, typename Sys_eqs>
class D<d,arg,Sys_eqs,CENTRAL>
{
	public:

87 88
	/*! \brief Calculate which colums of the Matrix are non zero
	 *
89 90
	 * \param g_map mapping grid
	 * \param kmap position where the derivative is calculated
91
	 * \param gs Grid info
92
	 * \param spacing grid spacing
93 94
	 * \param cols non-zero colums calculated by the function
	 * \param coeff coefficent (constant in front of the derivative)
incardon's avatar
incardon committed
95
	 *
96 97 98
	 * ### Example
	 *
	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
incardon's avatar
incardon committed
99 100
	 *
	 */
101 102 103 104 105 106
	inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map,
			                 grid_dist_key_dx<Sys_eqs::dims> & kmap,
							 const grid_sm<Sys_eqs::dims,void> & gs,
							 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
							 std::unordered_map<long int,typename Sys_eqs::stype > & cols,
							 typename Sys_eqs::stype coeff)
incardon's avatar
incardon committed
107
	{
108
		// if the system is staggered the CENTRAL derivative is equivalent to a forward derivative
109
		if (is_grid_staggered<Sys_eqs>::value())
110
		{
111
			D<d,arg,Sys_eqs,BACKWARD>::value(g_map,kmap,gs,spacing,cols,coeff);
112 113 114
			return;
		}

115 116
		long int old_val = kmap.getKeyRef().get(d);
		kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
117
		arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]/2.0 );
118
		kmap.getKeyRef().set_d(d,old_val);
incardon's avatar
incardon committed
119

120 121 122

		old_val = kmap.getKeyRef().get(d);
		kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
123
		arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]/2.0 );
124
		kmap.getKeyRef().set_d(d,old_val);
incardon's avatar
incardon committed
125 126 127 128 129
	}


	/*! \brief Calculate the position where the derivative is calculated
	 *
130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
	 * In case on non staggered case this function just return a null grid_key, in case of staggered,
	 *  it calculate how the operator shift in the cell
	 *
 	 	 \verbatim

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

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

		\endverbatim
	 *
	 * Consider this 2D staggered cell and a second order central derivative scheme, this lead to
incardon's avatar
incardon committed
148
	 *
149 150 151 152 153 154 155 156
	 * \f$ \frac{\partial v_y}{\partial x} \f$ is calculated on position (*), so the function return the grid_key {0,0}
	 *
	 * \f$ \frac{\partial v_y}{\partial y} \f$ is calculated on position (0), so the function return the grid_key {-1,-1}
	 *
	 * \f$ \frac{\partial v_x}{\partial x} \f$ is calculated on position (0), so the function return the grid_key {-1,-1}
	 *
	 * \f$ \frac{\partial v_x}{\partial y} \f$ is calculated on position (*), so the function return the grid_key {0,0}
	 *
157
	 * \param pos position where we are calculating the derivative
158 159
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
incardon's avatar
incardon committed
160
	 *
161 162
	 * \return where (in which cell grid) the derivative is calculated
	 *
incardon's avatar
incardon committed
163
	 */
164 165 166
	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos,
			                                          const grid_sm<Sys_eqs::dims,void> & gs,
													  const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
incardon's avatar
incardon committed
167
	{
168
		auto arg_pos = arg::position(pos,gs,s_pos);
incardon's avatar
incardon committed
169 170
		if (is_grid_staggered<Sys_eqs>::value())
		{
171
			if (arg_pos.get(d) == -1)
incardon's avatar
incardon committed
172
			{
173 174
				arg_pos.set_d(d,0);
				return arg_pos;
incardon's avatar
incardon committed
175 176 177
			}
			else
			{
178 179
				arg_pos.set_d(d,-1);
				return arg_pos;
incardon's avatar
incardon committed
180 181 182
			}
		}

183
		return arg_pos;
incardon's avatar
incardon committed
184 185 186 187
	}
};


188 189 190 191 192 193
/*! \brief Second order one sided Derivative scheme on direction i
 *
 * \verbatim
 *
 *  -1.5    2.0   -0.5
 *    +------*------*
incardon's avatar
incardon committed
194
 *
195 196 197 198 199 200 201 202 203 204 205
 * or
 *
 *  -0.5    2.0   -1.5
 *    *------*------+
 *
 *  in the bulk
 *
 *  -1        +1
 *   *---+---*
 *
 * \endverbatim
incardon's avatar
incardon committed
206 207 208 209 210 211 212
 *
 */
template<unsigned int d, typename arg, typename Sys_eqs>
class D<d,arg,Sys_eqs,CENTRAL_B_ONE_SIDE>
{
public:

213 214
	/*! \brief Calculate which colums of the Matrix are non zero
	 *
215 216
	 * \param g_map mapping grid points
	 * \param kmap position where the derivative is calculated
217
	 * \param gs Grid info
218
	 * \param spacing of the grid
219 220 221 222
	 * \param cols non-zero colums calculated by the function
	 * \param coeff coefficent (constant in front of the derivative)
	 *
	 * ### Example
incardon's avatar
incardon committed
223
	 *
224
	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
incardon's avatar
incardon committed
225 226
	 *
	 */
227 228 229 230 231 232
	static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map,
			          grid_dist_key_dx<Sys_eqs::dims> & kmap,
					  const grid_sm<Sys_eqs::dims,void> & gs,
					  typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
					  std::unordered_map<long int,typename Sys_eqs::stype > & cols,
					  typename Sys_eqs::stype coeff)
incardon's avatar
incardon committed
233 234 235
	{
#ifdef SE_CLASS1
		if (Sys_eqs::boundary[d] == PERIODIC)
236
			std::cerr << __FILE__ << ":" << __LINE__ << " error, it make no sense use one sided derivation with periodic boundary, please use CENTRAL\n";
incardon's avatar
incardon committed
237 238
#endif

239 240
		grid_key_dx<Sys_eqs::dims> pos = g_map.getGKey(kmap);

incardon's avatar
incardon committed
241 242
		if (pos.get(d) == (long int)gs.size(d)-1 )
		{
243
			arg::value(g_map,kmap,gs,spacing,cols,1.5*coeff/spacing[d]);
incardon's avatar
incardon committed
244

245 246
			long int old_val = kmap.getKeyRef().get(d);
			kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
247
			arg::value(g_map,kmap,gs,spacing,cols,-2.0*coeff/spacing[d]);
248
			kmap.getKeyRef().set_d(d,old_val);
incardon's avatar
incardon committed
249

250
			old_val = kmap.getKeyRef().get(d);
incardon's avatar
incardon committed
251
			kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 2);
252
			arg::value(g_map,kmap,gs,spacing,cols,0.5*coeff/spacing[d]);
253
			kmap.getKeyRef().set_d(d,old_val);
incardon's avatar
incardon committed
254 255 256
		}
		else if (pos.get(d) == 0)
		{
257
			arg::value(g_map,kmap,gs,spacing,cols,-1.5*coeff/spacing[d]);
incardon's avatar
incardon committed
258

259 260
			long int old_val = kmap.getKeyRef().get(d);
			kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
261
			arg::value(g_map,kmap,gs,spacing,cols,2.0*coeff/spacing[d]);
262
			kmap.getKeyRef().set_d(d,old_val);
incardon's avatar
incardon committed
263

264
			old_val = kmap.getKeyRef().get(d);
incardon's avatar
incardon committed
265
			kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 2);
266
			arg::value(g_map,kmap,gs,spacing,cols,-0.5*coeff/spacing[d]);
267
			kmap.getKeyRef().set_d(d,old_val);
incardon's avatar
incardon committed
268 269 270
		}
		else
		{
271 272
			long int old_val = kmap.getKeyRef().get(d);
			kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
273
			arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]);
274 275 276 277
			kmap.getKeyRef().set_d(d,old_val);

			old_val = kmap.getKeyRef().get(d);
			kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
278
			arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]);
279
			kmap.getKeyRef().set_d(d,old_val);
incardon's avatar
incardon committed
280 281 282 283 284
		}
	}

	/*! \brief Calculate the position where the derivative is calculated
	 *
285 286
	 * In case on non staggered case this function just return a null grid_key, in case of staggered,
	 *  it calculate how the operator shift in the cell
incardon's avatar
incardon committed
287
	 *
288
 	 	 \verbatim
incardon's avatar
incardon committed
289

290 291 292 293 294
		+--$--+
		|     |
		#  *  #
		|     |
		0--$--+
incardon's avatar
incardon committed
295

296 297 298 299 300 301 302 303 304 305
	  # = velocity(y)
	  $ = velocity(x)
	  * = pressure

		\endverbatim
	 *
	 * In the one side stencil the cell position depend if you are or not at the boundary.
	 * outside the boundary is simply the central scheme, at the boundary it is simply the
	 * staggered position of the property
	 *
306
	 * \param pos position where we are calculating the derivative
307 308 309
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
	 *
310 311
	 * \return where (in which cell grid) the derivative is calculated
	 *
312 313 314 315
	 */
	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
	{
		return arg::position(pos,gs,s_pos);
incardon's avatar
incardon committed
316 317 318 319
	}
};


320 321 322 323 324 325 326 327
/*! \brief First order FORWARD derivative on direction i
 *
 * \verbatim
 *
 *  -1.0    1.0
 *    +------*
 *
 * \endverbatim
328 329 330 331 332 333 334
 *
 */
template<unsigned int d, typename arg, typename Sys_eqs>
class D<d,arg,Sys_eqs,FORWARD>
{
	public:

335 336
	/*! \brief Calculate which colums of the Matrix are non zero
	 *
337 338
	 * \param g_map mapping grid
	 * \param kmap position where the derivative is calculated
339
	 * \param gs Grid info
340
	 * \param spacing grid spacing
341 342
	 * \param cols non-zero colums calculated by the function
	 * \param coeff coefficent (constant in front of the derivative)
343
	 *
344 345 346
	 * ### Example
	 *
	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
347 348
	 *
	 */
349 350 351 352 353 354
	inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map,
			                 grid_dist_key_dx<Sys_eqs::dims> & kmap,
							 const grid_sm<Sys_eqs::dims,void> & gs,
							 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
							 std::unordered_map<long int,typename Sys_eqs::stype > & cols,
							 typename Sys_eqs::stype coeff)
355
	{
356 357 358

		long int old_val = kmap.getKeyRef().get(d);
		kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
359
		arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]);
360
		kmap.getKeyRef().set_d(d,old_val);
361 362

		// backward
363
		arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]);
364 365 366 367 368
	}


	/*! \brief Calculate the position where the derivative is calculated
	 *
369 370
	 * In case of non staggered case this function just return a null grid_key, in case of staggered,
	 * the FORWARD scheme return the position of the staggered property
371
	 *
372
	 * \param pos position where we are calculating the derivative
373 374
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
375
	 *
376 377
	 * \return where (in which cell grid) the derivative is calculated
	 *
378
	 */
379 380 381
	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos,
			                                          const grid_sm<Sys_eqs::dims,void> & gs,
													  const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
382
	{
383
		return arg::position(pos,gs,s_pos);
384 385 386
	}
};

387 388 389 390 391 392 393 394
/*! \brief First order BACKWARD derivative on direction i
 *
 * \verbatim
 *
 *  -1.0    1.0
 *    *------+
 *
 * \endverbatim
395 396 397 398 399 400 401
 *
 */
template<unsigned int d, typename arg, typename Sys_eqs>
class D<d,arg,Sys_eqs,BACKWARD>
{
	public:

402 403
	/*! \brief Calculate which colums of the Matrix are non zero
	 *
404 405
	 * \param g_map mapping grid
	 * \param kmap position where the derivative is calculated
406
	 * \param gs Grid info
407
	 * \param spacing of the grid
408 409 410 411
	 * \param cols non-zero colums calculated by the function
	 * \param coeff coefficent (constant in front of the derivative)
	 *
	 * ### Example
412
	 *
413
	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
414 415
	 *
	 */
416 417 418 419 420 421
	inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map,
			                 grid_dist_key_dx<Sys_eqs::dims> & kmap,
							 const grid_sm<Sys_eqs::dims,void> & gs,
							 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
							 std::unordered_map<long int,typename Sys_eqs::stype > & cols,
							 typename Sys_eqs::stype coeff)
422 423 424 425
	{

		long int old_val = kmap.getKeyRef().get(d);
		kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
426
		arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]);
427 428 429
		kmap.getKeyRef().set_d(d,old_val);

		// forward
430
		arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]);
431 432 433 434 435
	}


	/*! \brief Calculate the position where the derivative is calculated
	 *
436 437
	 * In case of non staggered case this function just return a null grid_key, in case of staggered,
	 * the BACKWARD scheme return the position of the staggered property
438
	 *
439
	 * \param pos position where we are calculating the derivative
440 441
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
442
	 *
443 444
	 * \return where (in which cell grid) the derivative is calculated
	 *
445
	 */
446
	inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
447
	{
448
		return arg::position(pos,gs,s_pos);
449 450
	}
};
incardon's avatar
incardon committed
451 452

#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_ */