1/*
2 * Derivative.hpp
3 *
4 * Created on: Oct 5, 2015
5 * Author: Pietro Incardona
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_
9#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_
10
11
12#include "util/mathutil.hpp"
13#include "Vector/map_vector.hpp"
14#include "Grid/comb.hpp"
15#include "FiniteDifference/util/common.hpp"
16#include "util/util_num.hpp"
17#include <unordered_map>
18#include "FD_util_include.hpp"
19
20/*! \brief Derivative second order on h (spacing)
21 *
22 * \tparam d on which dimension derive
23 * \tparam Field which field derive
24 * \tparam impl which implementation
25 *
26 */
27template<unsigned int d, typename Field, typename Sys_eqs, unsigned int impl=CENTRAL>
28class D
29{
30 /*! \brief Calculate which colums of the Matrix are non zero
31 *
32 * \param pos position where the derivative is calculated
33 * \param gs Grid info
34 * \param cols non-zero colums calculated by the function
35 * \param coeff coefficent (constant in front of the derivative)
36 *
37 * ### Example
38 *
39 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
40 *
41 */
42 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)
43 {
44 std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
45 }
46
47 /*! \brief Calculate the position where the derivative is calculated
48 *
49 * In case of non staggered case this function just return a null grid_key, in case of staggered,
50 * it calculate how the operator shift the calculation in the cell
51 *
52 * \param pos position where we are calculating the derivative
53 * \param gs Grid info
54 * \param s_pos staggered position of the properties
55 *
56 * \return where (in which cell) the derivative is calculated
57 *
58 */
59 inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos,
60 const grid_sm<Sys_eqs::dims,void> & gs,
61 const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
62 {
63 std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
64
65 return pos;
66 }
67};
68
69/*! \brief Second order central Derivative scheme on direction i
70 *
71 * \verbatim
72 *
73 * -1 +1
74 * *---+---*
75 *
76 * \endverbatim
77 *
78 */
79template<unsigned int d, typename arg, typename Sys_eqs>
80class D<d,arg,Sys_eqs,CENTRAL>
81{
82 public:
83
84 /*! \brief Calculate which colums of the Matrix are non zero
85 *
86 * \param g_map mapping grid
87 * \param kmap position where the derivative is calculated
88 * \param gs Grid info
89 * \param spacing grid spacing
90 * \param cols non-zero colums calculated by the function
91 * \param coeff coefficent (constant in front of the derivative)
92 *
93 * ### Example
94 *
95 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
96 *
97 */
98 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,
99 grid_dist_key_dx<Sys_eqs::dims> & kmap,
100 const grid_sm<Sys_eqs::dims,void> & gs,
101 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
102 std::unordered_map<long int,typename Sys_eqs::stype > & cols,
103 typename Sys_eqs::stype coeff)
104 {
105 // if the system is staggered the CENTRAL derivative is equivalent to a forward derivative
106 if (is_grid_staggered<Sys_eqs>::value())
107 {
108 D<d,arg,Sys_eqs,BACKWARD>::value(g_map,kmap,gs,spacing,cols,coeff);
109 return;
110 }
111
112 long int old_val = kmap.getKeyRef().get(d);
113 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
114 arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]/2.0 );
115 kmap.getKeyRef().set_d(d,old_val);
116
117
118 old_val = kmap.getKeyRef().get(d);
119 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
120 arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]/2.0 );
121 kmap.getKeyRef().set_d(d,old_val);
122 }
123
124
125 /*! \brief Calculate the position where the derivative is calculated
126 *
127 * In case on non staggered case this function just return a null grid_key, in case of staggered,
128 * it calculate how the operator shift in the cell
129 *
130 \verbatim
131
132 +--$--+
133 | |
134 # * #
135 | |
136 0--$--+
137
138 # = velocity(y)
139 $ = velocity(x)
140 * = pressure
141
142 \endverbatim
143 *
144 * Consider this 2D staggered cell and a second order central derivative scheme, this lead to
145 *
146 * \f$ \frac{\partial v_y}{\partial x} \f$ is calculated on position (*), so the function return the grid_key {0,0}
147 *
148 * \f$ \frac{\partial v_y}{\partial y} \f$ is calculated on position (0), so the function return the grid_key {-1,-1}
149 *
150 * \f$ \frac{\partial v_x}{\partial x} \f$ is calculated on position (0), so the function return the grid_key {-1,-1}
151 *
152 * \f$ \frac{\partial v_x}{\partial y} \f$ is calculated on position (*), so the function return the grid_key {0,0}
153 *
154 * \param pos position where we are calculating the derivative
155 * \param gs Grid info
156 * \param s_pos staggered position of the properties
157 *
158 * \return where (in which cell grid) the derivative is calculated
159 *
160 */
161 inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos,
162 const grid_sm<Sys_eqs::dims,void> & gs,
163 const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
164 {
165 auto arg_pos = arg::position(pos,gs,s_pos);
166 if (is_grid_staggered<Sys_eqs>::value())
167 {
168 if (arg_pos.get(d) == -1)
169 {
170 arg_pos.set_d(d,0);
171 return arg_pos;
172 }
173 else
174 {
175 arg_pos.set_d(d,-1);
176 return arg_pos;
177 }
178 }
179
180 return arg_pos;
181 }
182};
183
184
185/*! \brief Second order one sided Derivative scheme on direction i
186 *
187 * \verbatim
188 *
189 * -1.5 2.0 -0.5
190 * +------*------*
191 *
192 * or
193 *
194 * -0.5 2.0 -1.5
195 * *------*------+
196 *
197 * in the bulk
198 *
199 * -1 +1
200 * *---+---*
201 *
202 * \endverbatim
203 *
204 */
205template<unsigned int d, typename arg, typename Sys_eqs>
206class D<d,arg,Sys_eqs,CENTRAL_B_ONE_SIDE>
207{
208public:
209
210 /*! \brief Calculate which colums of the Matrix are non zero
211 *
212 * \param g_map mapping grid points
213 * \param kmap position where the derivative is calculated
214 * \param gs Grid info
215 * \param spacing of the grid
216 * \param cols non-zero colums calculated by the function
217 * \param coeff coefficent (constant in front of the derivative)
218 *
219 * ### Example
220 *
221 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
222 *
223 */
224 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,
225 grid_dist_key_dx<Sys_eqs::dims> & kmap,
226 const grid_sm<Sys_eqs::dims,void> & gs,
227 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
228 std::unordered_map<long int,typename Sys_eqs::stype > & cols,
229 typename Sys_eqs::stype coeff)
230 {
231#ifdef SE_CLASS1
232 if (Sys_eqs::boundary[d] == PERIODIC)
233 std::cerr << __FILE__ << ":" << __LINE__ << " error, it make no sense use one sided derivation with periodic boundary, please use CENTRAL\n";
234#endif
235
236 grid_key_dx<Sys_eqs::dims> pos = g_map.getGKey(kmap);
237
238 if (pos.get(d) == (long int)gs.size(d)-1 )
239 {
240 arg::value(g_map,kmap,gs,spacing,cols,1.5*coeff/spacing[d]);
241
242 long int old_val = kmap.getKeyRef().get(d);
243 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
244 arg::value(g_map,kmap,gs,spacing,cols,-2.0*coeff/spacing[d]);
245 kmap.getKeyRef().set_d(d,old_val);
246
247 old_val = kmap.getKeyRef().get(d);
248 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 2);
249 arg::value(g_map,kmap,gs,spacing,cols,0.5*coeff/spacing[d]);
250 kmap.getKeyRef().set_d(d,old_val);
251 }
252 else if (pos.get(d) == 0)
253 {
254 arg::value(g_map,kmap,gs,spacing,cols,-1.5*coeff/spacing[d]);
255
256 long int old_val = kmap.getKeyRef().get(d);
257 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
258 arg::value(g_map,kmap,gs,spacing,cols,2.0*coeff/spacing[d]);
259 kmap.getKeyRef().set_d(d,old_val);
260
261 old_val = kmap.getKeyRef().get(d);
262 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 2);
263 arg::value(g_map,kmap,gs,spacing,cols,-0.5*coeff/spacing[d]);
264 kmap.getKeyRef().set_d(d,old_val);
265 }
266 else
267 {
268 long int old_val = kmap.getKeyRef().get(d);
269 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
270 arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]);
271 kmap.getKeyRef().set_d(d,old_val);
272
273 old_val = kmap.getKeyRef().get(d);
274 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
275 arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]);
276 kmap.getKeyRef().set_d(d,old_val);
277 }
278 }
279
280 /*! \brief Calculate the position where the derivative is calculated
281 *
282 * In case on non staggered case this function just return a null grid_key, in case of staggered,
283 * it calculate how the operator shift in the cell
284 *
285 \verbatim
286
287 +--$--+
288 | |
289 # * #
290 | |
291 0--$--+
292
293 # = velocity(y)
294 $ = velocity(x)
295 * = pressure
296
297 \endverbatim
298 *
299 * In the one side stencil the cell position depend if you are or not at the boundary.
300 * outside the boundary is simply the central scheme, at the boundary it is simply the
301 * staggered position of the property
302 *
303 * \param pos position where we are calculating the derivative
304 * \param gs Grid info
305 * \param s_pos staggered position of the properties
306 *
307 * \return where (in which cell grid) the derivative is calculated
308 *
309 */
310 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])
311 {
312 return arg::position(pos,gs,s_pos);
313 }
314};
315
316
317/*! \brief First order FORWARD derivative on direction i
318 *
319 * \verbatim
320 *
321 * -1.0 1.0
322 * +------*
323 *
324 * \endverbatim
325 *
326 */
327template<unsigned int d, typename arg, typename Sys_eqs>
328class D<d,arg,Sys_eqs,FORWARD>
329{
330 public:
331
332 /*! \brief Calculate which colums of the Matrix are non zero
333 *
334 * \param g_map mapping grid
335 * \param kmap position where the derivative is calculated
336 * \param gs Grid info
337 * \param spacing grid spacing
338 * \param cols non-zero colums calculated by the function
339 * \param coeff coefficent (constant in front of the derivative)
340 *
341 * ### Example
342 *
343 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
344 *
345 */
346 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,
347 grid_dist_key_dx<Sys_eqs::dims> & kmap,
348 const grid_sm<Sys_eqs::dims,void> & gs,
349 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
350 std::unordered_map<long int,typename Sys_eqs::stype > & cols,
351 typename Sys_eqs::stype coeff)
352 {
353
354 long int old_val = kmap.getKeyRef().get(d);
355 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
356 arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]);
357 kmap.getKeyRef().set_d(d,old_val);
358
359 // backward
360 arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]);
361 }
362
363
364 /*! \brief Calculate the position where the derivative is calculated
365 *
366 * In case of non staggered case this function just return a null grid_key, in case of staggered,
367 * the FORWARD scheme return the position of the staggered property
368 *
369 * \param pos position where we are calculating the derivative
370 * \param gs Grid info
371 * \param s_pos staggered position of the properties
372 *
373 * \return where (in which cell grid) the derivative is calculated
374 *
375 */
376 inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos,
377 const grid_sm<Sys_eqs::dims,void> & gs,
378 const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
379 {
380 return arg::position(pos,gs,s_pos);
381 }
382};
383
384/*! \brief First order BACKWARD derivative on direction i
385 *
386 * \verbatim
387 *
388 * -1.0 1.0
389 * *------+
390 *
391 * \endverbatim
392 *
393 */
394template<unsigned int d, typename arg, typename Sys_eqs>
395class D<d,arg,Sys_eqs,BACKWARD>
396{
397 public:
398
399 /*! \brief Calculate which colums of the Matrix are non zero
400 *
401 * \param g_map mapping grid
402 * \param kmap position where the derivative is calculated
403 * \param gs Grid info
404 * \param spacing of the grid
405 * \param cols non-zero colums calculated by the function
406 * \param coeff coefficent (constant in front of the derivative)
407 *
408 * ### Example
409 *
410 * \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
411 *
412 */
413 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,
414 grid_dist_key_dx<Sys_eqs::dims> & kmap,
415 const grid_sm<Sys_eqs::dims,void> & gs,
416 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
417 std::unordered_map<long int,typename Sys_eqs::stype > & cols,
418 typename Sys_eqs::stype coeff)
419 {
420
421 long int old_val = kmap.getKeyRef().get(d);
422 kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
423 arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]);
424 kmap.getKeyRef().set_d(d,old_val);
425
426 // forward
427 arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]);
428 }
429
430
431 /*! \brief Calculate the position where the derivative is calculated
432 *
433 * In case of non staggered case this function just return a null grid_key, in case of staggered,
434 * the BACKWARD scheme return the position of the staggered property
435 *
436 * \param pos position where we are calculating the derivative
437 * \param gs Grid info
438 * \param s_pos staggered position of the properties
439 *
440 * \return where (in which cell grid) the derivative is calculated
441 *
442 */
443 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])
444 {
445 return arg::position(pos,gs,s_pos);
446 }
447};
448
449#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_ */
450