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
#define CENTRAL_SYM 4
incardon's avatar
incardon committed
16 17 18 19 20

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

incardon's avatar
incardon committed
24 25 26 27 28 29 30 31 32 33
/*! \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
{
34 35 36 37 38 39
	/*! \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
40
	 *
41 42 43
	 * ### Example
	 *
	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
incardon's avatar
incardon committed
44 45
	 *
	 */
incardon's avatar
incardon committed
46
	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
47 48 49 50 51 52
	{
		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
	}

	/*! \brief Calculate the position where the derivative is calculated
	 *
53 54 55
	 * 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
	 *
56
	 * \param pos position where we are calculating the derivative
57 58
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
incardon's avatar
incardon committed
59
	 *
60 61
	 * \return where (in which cell) the derivative is calculated
	 *
incardon's avatar
incardon committed
62
	 */
63 64 65
	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
66 67
	{
		std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
68 69

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

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

88 89
	/*! \brief Calculate which colums of the Matrix are non zero
	 *
90 91
	 * \param g_map mapping grid
	 * \param kmap position where the derivative is calculated
92
	 * \param gs Grid info
93
	 * \param spacing grid spacing
94 95
	 * \param cols non-zero colums calculated by the function
	 * \param coeff coefficent (constant in front of the derivative)
incardon's avatar
incardon committed
96
	 *
97 98 99
	 * ### Example
	 *
	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
incardon's avatar
incardon committed
100 101
	 *
	 */
102 103 104 105 106 107
	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
108
	{
109
		// if the system is staggered the CENTRAL derivative is equivalent to a forward derivative
110
		if (is_grid_staggered<Sys_eqs>::value())
111
		{
112
			D<d,arg,Sys_eqs,BACKWARD>::value(g_map,kmap,gs,spacing,cols,coeff);
113 114 115
			return;
		}

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

121 122 123

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


	/*! \brief Calculate the position where the derivative is calculated
	 *
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148
	 * 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
149
	 *
150 151 152 153 154 155 156 157
	 * \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}
	 *
158
	 * \param pos position where we are calculating the derivative
159 160
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
incardon's avatar
incardon committed
161
	 *
162 163
	 * \return where (in which cell grid) the derivative is calculated
	 *
incardon's avatar
incardon committed
164
	 */
165 166 167
	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
168
	{
169
		auto arg_pos = arg::position(pos,gs,s_pos);
incardon's avatar
incardon committed
170 171
		if (is_grid_staggered<Sys_eqs>::value())
		{
172
			if (arg_pos.get(d) == -1)
incardon's avatar
incardon committed
173
			{
174 175
				arg_pos.set_d(d,0);
				return arg_pos;
incardon's avatar
incardon committed
176 177 178
			}
			else
			{
179 180
				arg_pos.set_d(d,-1);
				return arg_pos;
incardon's avatar
incardon committed
181 182 183
			}
		}

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


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

214 215
	/*! \brief Calculate which colums of the Matrix are non zero
	 *
216 217
	 * \param g_map mapping grid points
	 * \param kmap position where the derivative is calculated
218
	 * \param gs Grid info
219
	 * \param spacing of the grid
220 221 222 223
	 * \param cols non-zero colums calculated by the function
	 * \param coeff coefficent (constant in front of the derivative)
	 *
	 * ### Example
incardon's avatar
incardon committed
224
	 *
225
	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
incardon's avatar
incardon committed
226 227
	 *
	 */
228 229 230 231 232 233
	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
234 235 236
	{
#ifdef SE_CLASS1
		if (Sys_eqs::boundary[d] == PERIODIC)
237
			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
238 239
#endif

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

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

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

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

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

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

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

	/*! \brief Calculate the position where the derivative is calculated
	 *
286 287
	 * 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
288
	 *
289
 	 	 \verbatim
incardon's avatar
incardon committed
290

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

297 298 299 300 301 302 303 304 305 306
	  # = 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
	 *
307
	 * \param pos position where we are calculating the derivative
308 309 310
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
	 *
311 312
	 * \return where (in which cell grid) the derivative is calculated
	 *
313 314 315 316
	 */
	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
317 318 319 320
	}
};


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

336 337
	/*! \brief Calculate which colums of the Matrix are non zero
	 *
338 339
	 * \param g_map mapping grid
	 * \param kmap position where the derivative is calculated
340
	 * \param gs Grid info
341
	 * \param spacing grid spacing
342 343
	 * \param cols non-zero colums calculated by the function
	 * \param coeff coefficent (constant in front of the derivative)
344
	 *
345 346 347
	 * ### Example
	 *
	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
348 349
	 *
	 */
350 351 352 353 354 355
	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)
356
	{
357 358 359

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

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


	/*! \brief Calculate the position where the derivative is calculated
	 *
370 371
	 * 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
372
	 *
373
	 * \param pos position where we are calculating the derivative
374 375
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
376
	 *
377 378
	 * \return where (in which cell grid) the derivative is calculated
	 *
379
	 */
380 381 382
	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])
383
	{
384
		return arg::position(pos,gs,s_pos);
385 386 387
	}
};

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

403 404
	/*! \brief Calculate which colums of the Matrix are non zero
	 *
405 406
	 * \param g_map mapping grid
	 * \param kmap position where the derivative is calculated
407
	 * \param gs Grid info
408
	 * \param spacing of the grid
409 410 411 412
	 * \param cols non-zero colums calculated by the function
	 * \param coeff coefficent (constant in front of the derivative)
	 *
	 * ### Example
413
	 *
414
	 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
415 416
	 *
	 */
417 418 419 420 421 422
	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)
423 424 425 426
	{

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

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


	/*! \brief Calculate the position where the derivative is calculated
	 *
437 438
	 * 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
439
	 *
440
	 * \param pos position where we are calculating the derivative
441 442
	 * \param gs Grid info
	 * \param s_pos staggered position of the properties
443
	 *
444 445
	 * \return where (in which cell grid) the derivative is calculated
	 *
446
	 */
447
	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])
448
	{
449
		return arg::position(pos,gs,s_pos);
450 451
	}
};
incardon's avatar
incardon committed
452 453

#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_ */