OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
FD_op.hpp
1/*
2 * FD_op.hpp
3 *
4 * Created on: May 1, 2020
5 * Author: Abhinav Singh
6 */
7#ifndef FD_OP_HPP_
8#define FD_OP_HPP_
9#include "util/common.hpp"
10#include "FD_expressions.hpp"
11
12namespace FD
13{
14 constexpr int CENTRAL = 0;
15 constexpr int CENTRAL_ONE_SIDE_FORWARD = 1;
16 constexpr int CENTRAL_ONE_SIDE_BACKWARD = 2;
17 constexpr int CENTRAL_STAG = 3;
18 constexpr int CENTRAL_STAG_ONE_SIDE_FORWARD = 4;
19 constexpr int CENTRAL_STAG_ONE_SIDE_BACKWARD = 5;
20
21
22 // ord_d means order of the derivative and not order of convergence
23 template<unsigned int dir,unsigned int ord_d,unsigned int ord, unsigned int impl>
24 struct Derivative_impl
25 {
26 template<typename rtype,typename expr_type>
27 static rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key)
28 {
29 std::cerr << __FILE__ << ":" << __LINE__ << " Error the give derivative scheme with convergence order " << ord << " has not been implemented yet" << std::endl;
30 }
31 };
32
33
34 template<unsigned int dir>
35 struct Derivative_impl<dir,1,2,CENTRAL>
36 {
37 template<typename rtype,typename expr_type>
38 static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key,comb<expr_type::gtype::dims> & c_where)
39 {
40 // x0, dx are defined in proper dir є(x, y, z)
41 auto dx = o1.getGrid().spacing(dir);
42 long int x0 = key.getKeyRef().get(dir);
43
44 key.getKeyRef().set_d(dir, x0 + 1);
45 rtype ret = o1.value(key,c_where)*1.0/2.0;
46
47 key.getKeyRef().set_d(dir, x0 - 1);
48 ret += o1.value(key,c_where)*(-1.0/2.0);
49
50 key.getKeyRef().set_d(dir, x0);
51 return ret/dx;
52 }
53
54 template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
55 inline static void calculate_nz(o1_type o1,
56 const gmap_type & g_map,
57 grid_dist_key_dx<Sys_eqs::dims> & kmap,
58 const grid_sm<Sys_eqs::dims,void> & gs,
59 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
60 unordered_map_type & cols,
61 typename Sys_eqs::stype coeff,
62 unsigned int comp,
63 comb<Sys_eqs::dims> & c_where)
64 {
65 long int old_val = kmap.getKeyRef().get(dir);
66 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 1);
67 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,coeff/spacing[dir]/2.0,comp,c_where);
68 kmap.getKeyRef().set_d(dir,old_val);
69
70 old_val = kmap.getKeyRef().get(dir);
71 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 1);
72 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-coeff/spacing[dir]/2.0,comp,c_where);
73 kmap.getKeyRef().set_d(dir,old_val);
74 }
75 };
76
77 template<unsigned int dir>
78 struct Derivative_impl<dir,1,2,CENTRAL_STAG>
79 {
80 template<typename rtype,typename expr_type>
81 static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key,comb<expr_type::gtype::dims> & c_where)
82 {
83 rtype ret;
84 // x0, dx are defined in proper dir є(x, y, z)
85 auto dx = o1.getGrid().spacing(dir);
86 long int x0 = key.getKeyRef().get(dir);
87
88 if (c_where[dir] == -1)
89 {
90 char old_c = c_where[dir];
91 c_where[dir] = 0;
92 key.getKeyRef().set_d(dir, x0 - 1);
93 ret = o1.value(key,c_where)*(-1.0);
94
95 key.getKeyRef().set_d(dir, x0);
96 ret += o1.value(key,c_where)*(1.0);
97
98 c_where[dir] = old_c;
99 }
100 else
101 {
102 char old_c = c_where[dir];
103 c_where[dir] = -1;
104 key.getKeyRef().set_d(dir, x0);
105 ret = o1.value(key,c_where)*(-1.0);
106
107 key.getKeyRef().set_d(dir, x0 + 1);
108 ret += o1.value(key,c_where)*(1.0);
109
110 c_where[dir] = old_c;
111 }
112
113 return ret/dx;
114 }
115
116 template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
117 inline static void calculate_nz(o1_type o1,
118 const gmap_type & g_map,
119 grid_dist_key_dx<Sys_eqs::dims> & kmap,
120 const grid_sm<Sys_eqs::dims,void> & gs,
121 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
122 unordered_map_type & cols,
123 typename Sys_eqs::stype coeff,
124 unsigned int comp,
125 comb<Sys_eqs::dims> & c_where)
126 {
127 // x0, dx are defined in proper dir є(x, y, z)
128 long int old_val = kmap.getKeyRef().get(dir);
129
130 if (c_where[dir] == -1)
131 {
132 char old_c = c_where.c[dir];
133 c_where.c[dir] = 0;
134 kmap.getKeyRef().set_d(dir, old_val - 1);
135 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-coeff/spacing[dir],comp,c_where);
136
137 kmap.getKeyRef().set_d(dir, old_val);
138 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,coeff/spacing[dir],comp,c_where);
139
140 c_where.c[dir] = old_c;
141 }
142 else
143 {
144 char old_c = c_where.c[dir];
145 c_where.c[dir] = -1;
146 kmap.getKeyRef().set_d(dir, old_val);
147 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-coeff/spacing[dir],comp,c_where);
148
149 kmap.getKeyRef().set_d(dir, old_val + 1);
150 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,coeff/spacing[dir],comp,c_where);
151
152 c_where.c[dir] = old_c;
153 }
154
155 kmap.getKeyRef().set_d(dir,old_val);
156 }
157 };
158
159 template<unsigned int dir>
160 struct Derivative_impl<dir,2,2,CENTRAL>
161 {
162 template<typename rtype,typename expr_type>
163 static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
164 {
165 // x0, dx are defined in proper dir є(x, y, z)
166 auto dx = o1.getGrid().spacing(dir);
167 long int x0 = key.getKeyRef().get(dir);
168
169 rtype ret = o1.value(key,c_where)*(-2.0);
170
171 key.getKeyRef().set_d(dir, x0 + 1);
172 ret += o1.value(key,c_where)*(1.0);
173
174 key.getKeyRef().set_d(dir, x0 - 1);
175 ret += o1.value(key,c_where)*(1.0);
176
177 key.getKeyRef().set_d(dir, x0);
178 return ret/(dx*dx);
179 }
180
181 template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
182 inline static void calculate_nz(o1_type o1,
183 const gmap_type & g_map,
184 grid_dist_key_dx<Sys_eqs::dims> & kmap,
185 const grid_sm<Sys_eqs::dims,void> & gs,
186 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
187 unordered_map_type & cols,
188 typename Sys_eqs::stype coeff,
189 unsigned int comp,
190 comb<Sys_eqs::dims> & c_where)
191 {
192 long int old_val = kmap.getKeyRef().get(dir);
193 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-2.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
194
195 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 1);
196 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,coeff/(spacing[dir]*spacing[dir]),comp,c_where);
197 kmap.getKeyRef().set_d(dir,old_val);
198
199 old_val = kmap.getKeyRef().get(dir);
200 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 1);
201 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,coeff/(spacing[dir]*spacing[dir]),comp,c_where);
202 kmap.getKeyRef().set_d(dir,old_val);
203 }
204 };
205
206 enum one_side_direction
207 {
208 OS_CENTRAL,
209 OS_FORWARD,
210 OS_BACKWARD,
211 };
212
213 template<typename gtype, typename key_type>
214 one_side_direction use_one_side(gtype & g, unsigned int dir, key_type & key)
215 {
216 if (g.getDecomposition().periodicity()[dir] == true)
217 {
218 return one_side_direction::OS_CENTRAL;
219 }
220
221 auto keyg = g.getGKey(key);
222 if (keyg.get(dir) == 0)
223 {
224 return one_side_direction::OS_FORWARD;
225 }
226 else if (keyg.get(dir) == g.getGridInfoVoid().size(dir) - 1)
227 {
228 return one_side_direction::OS_BACKWARD;
229 }
230
231 return one_side_direction::OS_CENTRAL;
232 }
233
234 template<unsigned int dir>
235 struct Derivative_impl<dir,1,2,CENTRAL_ONE_SIDE_FORWARD>
236 {
237 template<typename rtype,typename expr_type>
238 static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
239 {
240 // x0, dx are defined in proper dir є(x, y, z)
241 auto dx = o1.getGrid().spacing(dir);
242 long int x0 = key.getKeyRef().get(dir);
243
244 rtype ret = o1.value(key,c_where)*(-3.0/2.0);
245
246 key.getKeyRef().set_d(dir, x0 + 1);
247 ret += o1.value(key,c_where)*2.0;
248
249 key.getKeyRef().set_d(dir, x0 + 2);
250 ret += o1.value(key,c_where)*(-0.5);
251
252 key.getKeyRef().set_d(dir, x0);
253 return ret/dx;
254 }
255
256 template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
257 inline static void calculate_nz(o1_type o1,
258 const gmap_type & g_map,
259 grid_dist_key_dx<Sys_eqs::dims> & kmap,
260 const grid_sm<Sys_eqs::dims,void> & gs,
261 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
262 unordered_map_type & cols,
263 typename Sys_eqs::stype coeff,
264 unsigned int comp,
265 comb<Sys_eqs::dims> & c_where)
266 { long int old_val = kmap.getKeyRef().get(dir);
267 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir));
268 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-3.0*coeff/spacing[dir]/2.0,comp,c_where);
269 kmap.getKeyRef().set_d(dir,old_val);
270 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 1);
271 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,2.0*coeff/spacing[dir],comp,c_where);
272 kmap.getKeyRef().set_d(dir,old_val);
273 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 2);
274 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-0.5*coeff/spacing[dir],comp,c_where);
275 kmap.getKeyRef().set_d(dir,old_val);
276 }
277 };
278
279 template<unsigned int dir>
280 struct Derivative_impl<dir,1,2,CENTRAL_STAG_ONE_SIDE_FORWARD>
281 {
282 template<typename rtype,typename expr_type>
283 static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
284 {
285
286 std::cout << __FILE__ << ":" << __LINE__ << " error we do not have implemented yet one sided staggered grid" << std::endl;
287
288 return 0.0;
289 }
290
291 template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
292 inline static void calculate_nz(o1_type o1,
293 const gmap_type & g_map,
294 grid_dist_key_dx<Sys_eqs::dims> & kmap,
295 const grid_sm<Sys_eqs::dims,void> & gs,
296 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
297 unordered_map_type & cols,
298 typename Sys_eqs::stype coeff,
299 unsigned int comp,
300 comb<Sys_eqs::dims> & c_where)
301 {
302 std::cout << __FILE__ << ":" << __LINE__ << " error we do not have implemented yet one sided staggered grid" << std::endl;
303 }
304 };
305
306 template<unsigned int dir>
307 struct Derivative_impl<dir,2,2,CENTRAL_ONE_SIDE_FORWARD>
308 {
309 template<typename rtype,typename expr_type>
310 static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
311 {
312 // x0, dx are defined in proper dir є(x, y, z)
313 auto dx = o1.getGrid().spacing(dir);
314 long int x0 = key.getKeyRef().get(dir);
315
316 rtype ret = o1.value(key,c_where)*(2.0);
317
318 key.getKeyRef().set_d(dir, x0 + 1);
319 ret += o1.value(key,c_where)*(-5.0);
320
321 key.getKeyRef().set_d(dir, x0 + 2);
322 ret += o1.value(key,c_where)*(4.0);
323
324 key.getKeyRef().set_d(dir, x0 + 3);
325 ret += o1.value(key,c_where)*(-1.0);
326
327 key.getKeyRef().set_d(dir, x0);
328 return ret/(dx*dx);
329 }
330
331 template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
332 inline static void calculate_nz(o1_type o1,
333 const gmap_type & g_map,
334 grid_dist_key_dx<Sys_eqs::dims> & kmap,
335 const grid_sm<Sys_eqs::dims,void> & gs,
336 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
337 unordered_map_type & cols,
338 typename Sys_eqs::stype coeff,
339 unsigned int comp,
340 comb<Sys_eqs::dims> & c_where)
341 {
342 long int old_val = kmap.getKeyRef().get(dir);
343 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir));
344 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,2.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
345 kmap.getKeyRef().set_d(dir,old_val);
346 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 1);
347 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-5.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
348 kmap.getKeyRef().set_d(dir,old_val);
349 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 2);
350 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,4.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
351 kmap.getKeyRef().set_d(dir,old_val);
352 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 3);
353 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-1.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
354 kmap.getKeyRef().set_d(dir,old_val);
355 }
356 };
357
358
359
360 template<unsigned int dir>
361 struct Derivative_impl<dir,1,2,CENTRAL_ONE_SIDE_BACKWARD>
362 {
363 template<typename rtype,typename expr_type>
364 static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
365 {
366 // x0, dx are defined in proper dir є(x, y, z)
367 auto dx = o1.getGrid().spacing(dir);
368 long int x0 = key.getKeyRef().get(dir);
369
370 rtype ret = o1.value(key,c_where)*(3.0/2.0);
371
372 key.getKeyRef().set_d(dir, x0 - 1);
373 ret += o1.value(key,c_where)*(-2.0);
374
375 key.getKeyRef().set_d(dir, x0 - 2);
376 ret += o1.value(key,c_where)*0.5;
377
378 key.getKeyRef().set_d(dir, x0);
379 return ret/dx;
380 }
381
382 template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
383 inline static void calculate_nz(o1_type o1,
384 const gmap_type & g_map,
385 grid_dist_key_dx<Sys_eqs::dims> & kmap,
386 const grid_sm<Sys_eqs::dims,void> & gs,
387 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
388 unordered_map_type & cols,
389 typename Sys_eqs::stype coeff,
390 unsigned int comp,
391 comb<Sys_eqs::dims> & c_where)
392 { long int old_val = kmap.getKeyRef().get(dir);
393 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir));
394 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,3.0*coeff/spacing[dir]/2.0,comp,c_where);
395 kmap.getKeyRef().set_d(dir,old_val);
396 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 1);
397 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-2.0*coeff/spacing[dir],comp,c_where);
398 kmap.getKeyRef().set_d(dir,old_val);
399 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 2);
400 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,0.5*coeff/spacing[dir],comp,c_where);
401 kmap.getKeyRef().set_d(dir,old_val);
402 }
403 };
404
405 template<unsigned int dir>
406 struct Derivative_impl<dir,1,2,CENTRAL_STAG_ONE_SIDE_BACKWARD>
407 {
408 template<typename rtype,typename expr_type>
409 static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
410 {
411
412 std::cout << __FILE__ << ":" << __LINE__ << " error we do not have implemented yet one sided staggered grid" << std::endl;
413
414 return 0.0;
415 }
416
417 template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
418 inline static void calculate_nz(o1_type o1,
419 const gmap_type & g_map,
420 grid_dist_key_dx<Sys_eqs::dims> & kmap,
421 const grid_sm<Sys_eqs::dims,void> & gs,
422 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
423 unordered_map_type & cols,
424 typename Sys_eqs::stype coeff,
425 unsigned int comp,
426 comb<Sys_eqs::dims> & c_where)
427 {
428 std::cout << __FILE__ << ":" << __LINE__ << " error we do not have implemented yet one sided staggered grid" << std::endl;
429 }
430 };
431
432 template<unsigned int dir>
433 struct Derivative_impl<dir,2,2,CENTRAL_ONE_SIDE_BACKWARD>
434 {
435 template<typename rtype,typename expr_type>
436 static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
437 {
438 // x0, dx are defined in proper dir є(x, y, z)
439 auto dx = o1.getGrid().spacing(dir);
440 long int x0 = key.getKeyRef().get(dir);
441
442 rtype ret = o1.value(key,c_where)*(2.0);
443
444 key.getKeyRef().set_d(dir, x0 - 1);
445 ret += o1.value(key,c_where)*(-5.0);
446
447 key.getKeyRef().set_d(dir, x0 - 2);
448 ret += o1.value(key,c_where)*(4.0);
449
450 key.getKeyRef().set_d(dir, x0 - 3);
451 ret += o1.value(key,c_where)*(-1.0);
452
453 key.getKeyRef().set_d(dir, x0);
454 return ret/(dx*dx);
455 }
456
457 template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
458 inline static void calculate_nz(o1_type o1,
459 const gmap_type & g_map,
460 grid_dist_key_dx<Sys_eqs::dims> & kmap,
461 const grid_sm<Sys_eqs::dims,void> & gs,
462 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
463 unordered_map_type & cols,
464 typename Sys_eqs::stype coeff,
465 unsigned int comp,
466 comb<Sys_eqs::dims> & c_where)
467 { long int old_val = kmap.getKeyRef().get(dir);
468 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir));
469 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,2.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
470 kmap.getKeyRef().set_d(dir,old_val);
471 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 1);
472 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-5.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
473 kmap.getKeyRef().set_d(dir,old_val);
474 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 2);
475 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,4.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
476 kmap.getKeyRef().set_d(dir,old_val);
477 kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 3);
478 o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-1.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
479 kmap.getKeyRef().set_d(dir,old_val);
480 }
481 };
482
483 template<unsigned dir_, unsigned int ord_d_, unsigned int ord_, unsigned int impl_>
484 struct GRID_DERIVATIVE
485 {
486 typedef boost::mpl::int_<dir_> dir;
487 typedef boost::mpl::int_<ord_> ord;
488 typedef boost::mpl::int_<ord_d_> ord_d;
489 typedef boost::mpl::int_<impl_> impl;
490 };
491
492 template <typename exp1, unsigned int dir, unsigned int ord_d,unsigned int ord, unsigned int impl>
493 class grid_dist_expression_op<exp1,void,GRID_DERIVATIVE<dir,ord_d,ord,impl> >
494 {
496 const exp1 o1;
497
498 public:
499
500 typedef typename exp1::gtype gtype;
501
503 inline grid_dist_expression_op(const exp1 & o1)
504 :o1(o1)
505 {}
506
512 inline void init() const
513 {
514 o1.init();
515 }
516
524 inline auto value(grid_dist_key_dx<gtype::dims> & key, comb<gtype::dims> & c_where) const -> typename std::remove_reference<decltype(o1.value(key,c_where))>::type
525 {
526 typedef typename std::remove_reference<decltype(o1.value(key,c_where))>::type r_type;
527
528 r_type val;
529
530 one_side_direction os = use_one_side(getGrid(),dir,key);
531
532 if (os == one_side_direction::OS_CENTRAL || getGrid().is_staggered() == true)
533 {
534 val = Derivative_impl<dir,ord_d,ord,impl>::template calculate<r_type>(o1,key,c_where);
535 }
536 else if (os == one_side_direction::OS_FORWARD)
537 {
538 val = Derivative_impl<dir,ord_d,ord,impl+1>::template calculate<r_type>(o1,key,c_where);
539 }
540 else
541 {
542 val = Derivative_impl<dir,ord_d,ord,impl+2>::template calculate<r_type>(o1,key,c_where);
543 }
544
545 return val;
546 }
547
548
549 template<typename Sys_eqs, typename gmap_type, typename unordered_map_type>
550 inline void value_nz(const gmap_type & g_map,
551 grid_dist_key_dx<Sys_eqs::dims> & kmap,
552 const grid_sm<Sys_eqs::dims,void> & gs,
553 typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
554 unordered_map_type & cols,
555 typename Sys_eqs::stype coeff,
556 unsigned int comp,
557 comb<Sys_eqs::dims> & c_where) const
558 {
559 //o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,coeff,comp);
560 long int old_val = kmap.getKeyRef().get(dir);
561 one_side_direction os = use_one_side(getGrid(),dir,kmap);
562 if (os == one_side_direction::OS_CENTRAL || getGrid().is_staggered() == true)
563 {
564 Derivative_impl<dir,ord_d,ord,impl>::template calculate_nz<Sys_eqs>(o1,g_map,kmap,gs,spacing,cols,coeff,comp,c_where);
565 }
566 else if (os == one_side_direction::OS_FORWARD)
567 {
568 Derivative_impl<dir,ord_d,ord,impl+1>::template calculate_nz<Sys_eqs>(o1,g_map,kmap,gs,spacing,cols,coeff,comp,c_where);
569 }
570 else
571 {
572 Derivative_impl<dir,ord_d,ord,impl+2>::template calculate_nz<Sys_eqs>(o1,g_map,kmap,gs,spacing,cols,coeff,comp,c_where);
573 }
574
575 // cols[g_map.template getProp<0>(kmap)*Sys_eqs::nvar + FD::var_id + comp] += coeff;
576 }
577
578
586 gtype & getGrid()
587 {
588 return o1.getGrid();
589 }
590
598 const gtype & getGrid() const
599 {
600 return o1.getGrid();
601 }
602 };
603
604
605
606
607 template<unsigned int dir,unsigned int ord_d ,unsigned int ord, unsigned int impl>
608 class Derivative
609 {
610
611 public:
612
613 Derivative()
614 {
615 }
616
617 template<typename operand_type>
618 grid_dist_expression_op<operand_type,void,GRID_DERIVATIVE<dir,ord_d,ord,impl>> operator()(operand_type arg)
619 {
620 return grid_dist_expression_op<operand_type,void,GRID_DERIVATIVE<dir,ord_d,ord,impl>>(arg);
621 }
622 };
623
624 template<unsigned int dim, unsigned int impl>
625 class Laplacian
626 {
627
628 public:
629
630 Laplacian()
631 {
632 }
633
634 template<typename operand_type>
635 FD::grid_dist_expression_op<
636 FD::grid_dist_expression_op<operand_type,void,GRID_DERIVATIVE<0,2,2,impl>>,
637 FD::grid_dist_expression_op<operand_type,void,GRID_DERIVATIVE<1,2,2,impl>>,
638 FD::sum>
639 operator()(operand_type arg)
640 {
641 grid_dist_expression_op<operand_type,void,GRID_DERIVATIVE<0,2,2,impl>> d_xx(arg);
642 grid_dist_expression_op<operand_type,void,GRID_DERIVATIVE<1,2,2,impl>> d_yy(arg);
643 return d_xx + d_yy;
644 }
645 };
646
647 class L2Error
648 {
649 public:
650 L2Error() {}
651
652 template<unsigned int p1, unsigned int p2, typename g1, typename g2, unsigned int impl_p1, unsigned int impl_p2>
653 auto operator()(const FD::grid_dist_expression<p1,g1,impl_p1> ga, const FD::grid_dist_expression<p2,g2,impl_p2> gb) const ->
654 typename std::remove_reference<decltype(ga.getGrid().template getProp<p1>(ga.getGrid().getDomainIterator().get()))>::type
655 {
656 typename std::remove_reference<decltype(ga.getGrid().template getProp<p1>(ga.getGrid().getDomainIterator().get()))>::type result = 0.0;
657 comb<FD::grid_dist_expression<p1,g1,impl_p1>::gtype::dims> s_pos;
658 auto diff_expr = ga + -1*gb;
659
660 auto & v_cl = create_vcluster();
661
662 auto it = ga.getGrid().getDomainIterator();
663 while (it.isNext())
664 {
665 auto key = it.get();
666 auto diff = diff_expr.value(key,s_pos);
667 result += diff*diff;
668 ++it;
669 }
670
671 v_cl.sum(result);
672 v_cl.execute();
673
674 return result;
675 }
676 };
677
678
679 class LInfError
680 {
681 public:
682 LInfError() {}
683
684 template<unsigned int p1, unsigned int p2, typename g1, typename g2, unsigned int impl_p1, unsigned int impl_p2>
685 auto operator()(const FD::grid_dist_expression<p1,g1,impl_p1> ga, const FD::grid_dist_expression<p2,g2,impl_p2> gb) const ->
686 typename std::remove_reference<decltype(ga.getGrid().template getProp<p1>(ga.getGrid().getDomainIterator().get()))>::type
687 {
688 typename std::remove_reference<decltype(ga.getGrid().template getProp<p1>(ga.getGrid().getDomainIterator().get()))>::type result = 0.0;
689 comb<FD::grid_dist_expression<p1,g1,impl_p1>::gtype::dims> s_pos;
690 auto diff_expr = ga + -1*gb;
691
692 auto & v_cl = create_vcluster();
693
694 auto it = ga.getGrid().getDomainIterator();
695 while (it.isNext())
696 {
697 auto key = it.get();
698 auto diff = fabs(diff_expr.value(key,s_pos));
699
700 if (diff > result) result = diff;
701 ++it;
702 }
703
704 v_cl.max(result);
705 v_cl.execute();
706
707 return result;
708 }
709 };
710
711 typedef Derivative<0,1,2,CENTRAL_STAG> Derivative_x_stag;
712 typedef Derivative<1,1,2,CENTRAL_STAG> Derivative_y_stag;
713 typedef Derivative<0,1,2,CENTRAL> Derivative_x;
714 typedef Derivative<1,1,2,CENTRAL> Derivative_y;
715 typedef Derivative<0,2,2,CENTRAL> Derivative_xx;
716 typedef Derivative<1,2,2,CENTRAL> Derivative_yy;
717 typedef Laplacian<2,CENTRAL> Lap;
718};
719
720
721#endif
FD::grid_dist_expression_op< exp1, void, GRID_DERIVATIVE< dir, ord_d, ord, impl > >::getGrid
gtype & getGrid()
Return the vector on which is acting.
Definition FD_op.hpp:586
FD::grid_dist_expression_op< exp1, void, GRID_DERIVATIVE< dir, ord_d, ord, impl > >::init
void init() const
This function must be called before value.
Definition FD_op.hpp:512
FD::grid_dist_expression_op< exp1, void, GRID_DERIVATIVE< dir, ord_d, ord, impl > >::grid_dist_expression_op
grid_dist_expression_op(const exp1 &o1)
Costruct a FD expression out of two expressions.
Definition FD_op.hpp:503
FD::grid_dist_expression_op< exp1, void, GRID_DERIVATIVE< dir, ord_d, ord, impl > >::getGrid
const gtype & getGrid() const
Return the vector on which is acting.
Definition FD_op.hpp:598
FD::grid_dist_expression_op< exp1, void, GRID_DERIVATIVE< dir, ord_d, ord, impl > >::o1
const exp1 o1
expression 1
Definition FD_op.hpp:496
FD::grid_dist_expression_op< exp1, void, GRID_DERIVATIVE< dir, ord_d, ord, impl > >::value
auto value(grid_dist_key_dx< gtype::dims > &key, comb< gtype::dims > &c_where) const -> typename std::remove_reference< decltype(o1.value(key, c_where))>::type
Evaluate the expression.
Definition FD_op.hpp:524
FD::grid_dist_expression_op
Definition FD_expressions.hpp:45
FD::grid_dist_expression
Definition FD_expressions.hpp:508
grid_dist_key_dx
Grid key for a distributed grid.
Definition grid_dist_key.hpp:44
grid_dist_key_dx::getKeyRef
base_key & getKeyRef()
Get the reference key.
Definition grid_dist_key.hpp:91
grid_sm
Declaration grid_sm.
Definition grid_sm.hpp:167
FD::Derivative_impl
Definition FD_op.hpp:25
FD::GRID_DERIVATIVE
Definition FD_op.hpp:485
comb
Position of the element of dimension d in the hyper-cube of dimension dim.
Definition comb.hpp:35
comb::c
signed char c[dim]
Array that store the combination.
Definition comb.hpp:37